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12 Water Quality Modelling and Prediction 


The most fundamental human needs for water are for drinking, cooking and 
personal hygiene. To meet these needs, the quality of the water used must pose no 
risk to human health. The quality of the water in nature also affects the condition 
of ecosystems that all living organisms depend on. At the same time, humans use 
water bodies as convenient sinks for the disposal of domestic, industrial and 
agricultural wastewaters. This of course degrades the quality of those water 
bodies. Water resources management involves the monitoring and management of 
water quality as much as the monitoring and management of water quantity. 
Various models have been developed to assist in predicting the water quality 
impacts of alternative land and water management policies and practices. This 
chapter introduces some of the main principles of water quality modelling. 


1. Introduction 

Water quality management is a critical component of 
overall integrated water resources management. Most 
users of water depend on adequate levels of water quality. 
When these levels are not met, these water users must 
either pay an additional cost for water treatment or incur 
at least increased risks of damage or loss. As populations 
and economies grow, more pollutants are generated. Many 
of these are waterborne, and hence can end up in surface 
and groundwater bodies. Increasingly, the major efforts 
and costs involved in water management are devoted to 
water quality protection and management. Conflicts 
among various users of water are increasingly over issues 
involving water quality as well as water quantity. 

Natural water bodies are able to serve many uses, 
including the transport and assimilation of waterborne 
wastes. But as natural water bodies assimilate these 
wastes, their quality changes. If the quality drops to 
the extent that other beneficial uses are adversely affected, 
the assimilative capacities of those water bodies have been 
exceeded with respect to those affected uses. Water 
quality management measures are actions taken to ensure 


that the total pollutant loads discharged into receiving 
water bodies do not exceed the ability of those water 
bodies to assimilate those loads while maintaining the 
levels of quality specified by quality standards set for 
those waters. 

What uses depend on water quality? One can identify 
almost any use. All living organisms require water of 
sufficient quantity and quality to survive, although differ- 
ent aquatic species can tolerate different levels of water 
quality. Regrettably, in most parts of the developed world 
it is no longer safe to drink natural surface or ground 
waters; they usually need to be treated before they become 
fit for human consumption. Treatment is not a practical 
option for recreational bathing, or for maintaining the 
health of fish, shellfish and other organisms found in 
natural aquatic ecosystems. Thus, standards specifying 
minimum acceptable levels of quality are set for most 
ambient waters. Various other uses have their own stan- 
dards as well. Irrigation water must not be too saline or 
contain toxic substances that can be absorbed by the 
plants or destroy microorganisms in the soil. Water qual- 
ity standards for industry can be very demanding, depend- 
ing of course on the particular industrial processes. 
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Pollutant loadings degrade water quality. High 
domestic wasteloads can result in high concentrations of 
bacteria, viruses and other organisms that affect human 
health. High organic loadings can reduce dissolved 
oxygen to levels that are fatal to parts of the aquatic 
ecosystem and cause obnoxious odours. Nutrient 
loadings from both urban and agricultural land runoff can 
cause excessive algae growth, which in turn may degrade 
the water aesthetically and recreationally, and upon 
death result in low dissolved oxygen levels. Toxic heavy 
metals and other micro-pollutants can accumulate in 
the bodies of aquatic organisms, including fish, making 
them unfit for human consumption even if they them- 
selves survive. 

Pollutant discharges originate from point and non- 
point sources. A common approach to controlling point 
source discharges, such as those from stormwater outfalls, 
municipal wastewater treatment plants or industries, is to 
impose standards specifying maximum allowable pollu- 
tant loads or concentrations in their effluents. This is 
often done in ways that are not economically efficient 
or even environmentally effective. Effluent standards 
typically do not take into account the particular assimila- 
tive capacities of the receiving water body. 

Non-point sources such as agricultural runoff or atmos- 
pheric deposition are less easily controlled, and hence it is 
difficult to apply effluent standards to non-point source 
pollutants, and their loadings can be much more signifi- 
cant than point source loadings. Management of non-point 
water quality impacts requires a more ambient-focused 
water quality management programme. 

The goal of an ambient water quality management 
programme is to establish appropriate standards for 
water quality in water bodies receiving pollutant loads, 
and then to ensure that these standards are met. Realistic 
standard-setting takes into account the basin’s hydro- 
logical, ecological and land use conditions, the potential 
uses of the receiving water body, and the institutional 
capacity to set and enforce water quality standards. 

Ambient-based water quality prediction and manage- 
ment involves considerable uncertainty. No one can 
predict what pollutant loadings will occur in the future, 
especially from area-wide non-point sources. In addition 
to uncertainties inherent in measuring the attainment 
of water quality standards, there are uncertainties in 
models used to determine sources of pollution, to allocate 


pollutant loads, and to predict the effectiveness of actions 
taken to meet water quality standards. The models 
available to help managers predict water quality impacts 
are relatively simple compared with the complexities of 
actual water systems. These limitations and uncertainties 
should be understood and addressed as water quality 
management decisions are based on their outputs. 

2. Establishing Ambient Water 
Quality Standards 

Identifying the intended uses of a water body - whether 
a lake, a section of a stream or an estuary - is a first step 
in setting water quality standards for that body. The most 
restrictive of the specific desired uses of a water body is 
termed a designated use. Barriers to achieving the desig- 
nated use are the presence of pollutants, or hydrological 
and geomorphic changes that affect the water quality. 

The designated use dictates the appropriate type 
of water quality standard. For example, a designated use 
of human recreation should protect humans from 
exposure to microbial pathogens while swimming, 
wading or boating. Other uses include those designed to 
protect humans and wildlife from consuming harmful 
substances in water, in fish and in shellfish. Aquatic-life 
uses include the protection and propagation of fish, 
shellfish and wildlife resources. 

Standards set upstream may affect the uses of water 
downstream. For example, small headwater streams 
may have aesthetic value but may not be able to support 
extensive recreational uses. However, their condition 
may affect the ability of a downstream area to achieve a 
particular designated use such as Tishable’ or ‘swim- 
mable’. In this case, the designated use for the smaller 
upstream water body may be defined in terms of 
achieving the designated use of the larger downstream 
water body. 

In many areas, human activities have altered the land- 
scape and aquatic ecosystems to the point where they 
cannot be restored to their pre-disturbance condition. For 
example, someone’s desire to establish a trout fish-farm in 
downtown Paris, Phnom Penh, Prague or Pretoria may 
not be attainable because of the development history of 
these areas or the altered hydrological regimes of the rivers 
flowing through them. Similarly, someone might wish to 
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designate an area near the outfall of a sewage treatment 
plant for shellfish harvesting, but health considerations 
would preclude any such use. Ambient water quality 
standards must be realistic. 

Designating the appropriate use for a water body is a 
policy decision that can be informed by the use of water 
quality prediction models of the type discussed in this 
chapter. However, the final standard selection should 
reflect a social consensus made while bearing in mind the 
current condition of the watershed, its pre-disturbance 
condition, the advantages derived from a certain desig- 
nated use, and the costs of achieving that use. 

2.1. Water-Use Criteria 

The designated use is a qualitative description of 
the desired condition of a water body. A criterion is a 
measurable indicator surrogate for use attainment. The 
criterion may be positioned at any point in the causal 
chain of boxes shown in Figure 12.1. 

In Box 1 of Figure 12.1 are measures of the pollutant 
discharge from a treatment plant (such as biological 
oxygen demand, ammonia (NH 3 ), pathogens and sus- 
pended sediments) or the amount of a pollutant entering 
the edge of a stream from runoff. A criterion at this posi- 
tion is referred to as an effluent standard. Criteria in Boxes 
2 and 3 are possible measures of ambient water quality 
conditions. Box 2 includes measures of a water quality 
parameter such as dissolved oxygen (DO), pH, total phos- 
phorus concentration, suspended sediment or tempera- 
ture. Criteria closer to the designated use (e.g. Box 3) 



Figure 12.1. Factors considered when determining designated 
use and associated water quality standards. 


include more combined or comprehensive measures 
of the biological community as a whole, such as the 
condition of the algal community (chlorophyll a) or 
a measure of contaminant concentration in fish tissue. 
Box 4 represents criteria that are associated with sources 
of pollution other than pollutants. These criteria might 
include measures such as flow timing and pattern 
(a hydrological criterion), abundance of non-indigenous 
taxa, or some quantification of channel modification 
(e.g. a decrease in sinuosity) (NRC, 2001). 

The more precise the statement of the designated use, 
the more accurate the criterion will be as an indicator of 
that use. For example, the criterion of fecal coliform count 
may be a suitable criterion for water contact recreation. 
The maximum allowable count itself may differ among 
water bodies that have water contact as their designated 
use, however. 

Surrogate indicators are often selected for use as 
criteria because they are easy to measure and in some 
cases are politically appealing. Although a surrogate 
indicator may have these appealing attributes, its useful- 
ness can be limited unless it can be logically related to 
a designated use. 

As with setting designated uses, the connections among 
water bodies and segments must be considered when 
determining criteria. For example, where a segment of a 
water body is designated as a mixing zone for a pollutant 
discharge, the criterion adopted should assure that the 
mixing zone use will not adversely affect the surrounding 
water body uses. Similarly, the desired condition of a small 
headwater stream may need to be specified in relation to 
other water bodies downstream; thus, an ambient nutrient 
criterion may be set in a small headwater stream to 
ensure a designated use in a downstream estuary, even if 
there are no local adverse impacts resulting from the nutri- 
ents in the small headwater stream, as previously 
discussed. Conversely, a high fecal coliform criterion 
may be permitted upstream of a recreational area if the 
fecal load dissipates before the flow reaches that area. 

3. Water Quality Model Use 

Monitoring data are the preferred form of information 
for identifying impaired waters (Appendix B). Model 
predictions might be used in addition to or instead of 
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monitoring data for several reasons: 

• Modelling might be feasible in some situations where 
monitoring is not. 

• Integrated monitoring and modelling systems could 
provide better information than one or the other 
alone for the same total cost. For example, regression 
analyses that correlate pollutant concentration with 
some more easily measurable factor (such as stream- 
flow) could be used to extend monitoring data for 
preliminary listing (of impared status) purposes. 
Models can also be used in a Bayesian framework 
to determine preliminary probability distributions of 
impairment that can help direct monitoring efforts and 
reduce the quantity of monitoring data needed for 
making listing decisions at a given level of reliability 
(see Chapter 7). 

• Modelling can be used to assess (predict) future water 
quality situations resulting from different management 
strategies. For example, assessing the improvement in 
water quality after a new wastewater treatment plant is 
built, or the effect of increased industrial growth and 
effluent discharges. 

A simple but useful modelling approach that may be 
used in the absence of monitoring data is ‘dilution calcula- 
tions’. In this approach the rate of pollutant loading from 
point sources in a water body is divided by the streamflow 
to give a set of estimated pollutant concentrations that may 
be compared to the standard. Simple dilution calculations 
assume conservative movement of pollutants. Thus, the 
use of dilution calculations will tend to be conservative 
and predict higher than actual concentrations for decaying 
pollutants. Of course, one could include a best estimate of 
the effects of decay processes in the dilution model. 

Combined runoff and water quality prediction models 
link stressors (sources of pollutants and pollution) to 
responses. Stressors include human activities likely to 
cause impairment, such as the presence of impervious 
surfaces in a watershed, cultivation of fields close to the 
stream, over-irrigation of crops with resulting polluted 
return flows, the discharge of domestic and industrial 
effluents into water bodies, installing dams and other 
channelization works, introduction of non-indigenous 
taxa and over-harvesting of fish. Indirect effects of humans 
include land cover changes that alter the rates of 
delivery of water, pollutants and sediment to water bodies. 


A review of direct and indirect effects of human 
activities suggests five major types of environmental 
stressors: 

• alterations in physical habitat 

• modifications in the seasonal flow of water 

• changes in the food base of the system 

• changes in interactions within the stream biota 

• release of contaminants (conventional pollutants) 
(Karr, 1990; NRC, 1992, 2001). 

Ideally, models designed to manage water quality should 
consider all five types of alternative management 
measures. A broad-based approach that considers these 
five features provides a more integrative approach to 
reduce the cause or causes of degradation (NRC, 1992). 

Models that relate stressors to responses can be of 
varying levels of complexity. Sometimes, they are sim- 
ple qualitative conceptual representations of the 
relationships among important variables and indicators 
of those variables, such as the statement ‘human activi- 
ties in a watershed affect water quality, including the 
condition of the river biota’. More quantitative models 
can be used to make predictions about the assimilative 
capacity of a water body, the movement of a pollutant 
from various point and non-point sources through a 
watershed, or the effectiveness of certain best manage- 
ment practices. 

3.1. Model Selection Criteria 

Water quality predictive models include both mathe- 
matical expressions and expert scientific judgement. They 
include process-based (mechanistic) models and 
data-based (statistical) models. The models should link 
management options to meaningful response variables 
(such as pollutant sources and water quality standard 
parameters). They should incorporate the entire ‘chain’ 
from stressors to responses. Process-based models should 
be consistent with scientific theory. Model prediction 
uncertainty should be reported. This provides decision- 
makers with estimates of the risks of options. To do this 
requires prediction error estimates (Chapter 9). 

Water quality management models should be appro- 
priate to the complexity of the situation and to the 
available data. Simple water quality problems can be 
addressed with simple models, while complex ones may 
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or may not require the use of more complex models. 
Models requiring large amounts of monitoring data 
should not be used in situations where such data are 
unavailable. Models should be flexible enough to allow 
updates and improvements as appropriate based on new 
research and monitoring data. 

Stakeholders need to accept the models proposed 
for use in any water quality management study. Given 
the increasing role of stakeholders in water management 
decision processes, they need to understand and accept 
the models being used, at least to the extent they wish to 
do so. Finally, the cost of maintaining and updating the 
model over time must be acceptable. 

Although predictions are typically made with the aid 
of mathematical models, there are certainly situations 
where expert judgement can be just as good. Reliance 
on professional judgement and simpler models is often 
acceptable, especially when data are limited. 

Highly detailed models require more time and 
are more expensive to develop and apply. Effective and 
efficient modelling for water quality management may 
dictate the use of simpler models. Complex modelling 
studies should be undertaken only if warranted by the 
complexity of the management problem. More complex 
modelling will not necessarily ensure that uncertainty is 
reduced, and in fact added complexity can compound 
problems of uncertainty analyses (Chapter 9). 

Placing a priority on process description usually leads 
to the development and use of complex mechanistic 
models rather than simpler mechanistic or empirical 
models. In some cases this may result in unnecessarily 
costly analyses. In addition, physical, chemical and bio- 
logical processes in terrestrial and aquatic environments 
are far too complex to be fully represented in even 
the most complicated models. For water quality manage- 
ment, the primary purpose of modelling should be to 
support decision-making. The inability to describe all 
relevant processes completely contributes to the uncer- 
tainty in the model predictions. 

3.2. Model Chains 

Many water quality management analyses require the use 
of a sequence of models, one feeding data into another. 
For example, consider the sequence or chain of models 
required for the prediction of fish and shellfish survival as 


a function of nutrient loadings into an estuary. Of interest 
to the stakeholders are the conditions of the fish and 
shellfish. One way to maintain healthy fish and shellfish 
stocks is to maintain sufficient levels of oxygen in the 
estuary. The way to do this is to control algae blooms. 
This in turn requires limits on the nutrient loadings to the 
estuary that can cause algae blooms, and subsequent 
dissolved oxygen deficits. The modelling challenge is to 
link nutrient loading to fish and shellfish survival. 

The negative effects of excessive nutrients (e.g. nitro- 
gen) in an estuary are shown in Figure 12.2. Nutrients 
stimulate the growth of algae. Algae die and accumulate 
on the bottom, where bacteria consume them. Under 
calm wind conditions density stratification occurs. 
Oxygen is depleted in the lower levels of water. Fish 
and shellfish may die or become weakened and more 
vulnerable to disease. 

A model consisting of a sequence of conditional 
probabilities can be defined to predict the probability of 
shellfish and fish abundance based on upstream nutrient 
loadings into the estuary that might cause problems 
for fish and shellfish populations. These conditional 
probabilities can be judgemental, mechanistic and/or sta- 
tistical. Each conditional probability can be a separate 
sub-model. Assuming each sub-model can identify a 
conditional probability distribution, the probability Pr{C 
I N} of a specified amount of carbon, C, given some spec- 
ified loading of a nutrient, say nitrogen, N, equals the 
probability Pr{C I A} of that given amount of carbon, 
given a concentration of algae biomass, A, times the prob- 
ability Pr{A | N, R} of that concentration of algae biomass 
given the nitrogen loading, N, and the river flow, R, 
times the probability Pr{R} of the river flow, R. In other 
words: 

Pr{C | N} = Pr{C I A}Pr{A | N, R}Pr{R} (12.1) 

An empirical process-based model of the type to be 
presented later in this chapter could be used to predict 
the concentration of algae and the chlorophyll violations 
on the basis of the river flow and nitrogen loadings. It 
could similarly predict the production of carbon, on the 
basis of algae biomass. A seasonal statistical regression 
model might be used to predict the likelihood of algae 
blooms based on algal biomass. A cross-system compari- 
son may be made to predict sediment oxygen demand. 
A relatively simple hydraulic model could be used to 
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Figure 12.2. The negative 
impacts of excessive nutrients 
in an estuary (NRC P 2001). 
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predict the duration of stratification and the frequency 
of hypoxia, given both the stratification duration and 
sediment oxygen demand. Expert judgement and fish 
survival models could be used to predict the shellfish 
abundance and fishkill and fish health probabilities. 

The biological endpoints ‘shell-fish survival’ and 
‘number of fishkills’, are meaningful indicators to stake- 
holders and can easily be related to designated water body 
use. Models and even conditional probabilities assigned 
to each link of the network in Figure 12.3 can reflect 
a combination of simple mechanisms, statistical (regres- 
sion) fitting and expert judgement. 

Advances in the mechanistic modelling of aquatic 
ecosystems have enabled us to include greater process 
(especially trophic) detail and complexity, as well as to 
perform dynamic simulations, although mechanistic 
ecosystem models have not advanced to the point of 
being able to predict community structure or biotic 
integrity. In this chapter, only some of the simpler 
mechanistic models will be introduced. More detail can 
be found in books solely devoted to water quality 


modelling (Chapra, 1997; McCutcheon, 1989; Orlob, 
1983; Schnoor, 1996; Thomann and Mueller, 1987) as 
well as the current professional journal literature. 

3.3. Model Data 

Data availability and accuracy are sources of concern 
in the development and use of models for water quality 
management. The complexity of models used for water 
quality management should be compatible with the 
quantity and quality of available data. The use of complex 
mechanistic models for water quality prediction in 
situations with little useful water quality data does not 
compensate for that lack of data. Model complexity 
can give the impression of credibility, but this is usually 
misleading. 

It is often preferable to begin with simple models and 
then, over time, add additional complexity as justified 
by the collection and analysis of additional data. This 
strategy makes efficient use of resources. It targets the 
effort toward information and models that will reduce 
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Figure 12.3. Cause and effect diagram for 
estuary eutrophication due to excessive 
nutrient loadings (Borsuk et al. p 2004). 


the uncertainty as the analysis proceeds. Models should 
be selected (simple versus complex) in part on the basis 
of the data available to support their use. 

Water quality models of water bodies receiving 
pollutant discharges require those pollutant loadings as 
input data. These pollutant discharges can be from point 
and non-point sources. Point source discharges are much 
easier to measure, monitor and estimate than non-point 
source inputs. Non-point discharge data often come 
from rainfall-runoff models that attempt to predict the 
quantity of runoff and its constituent concentrations. The 
reliability of the predictions from these models is not very 
good, especially if short time periods (e.g. each day or 
week) are being simulated. Their average values over 
longer time periods (e.g. a month or year) tend to be 
more reliable. This is mainly because the short-term 


inputs to those models, such as constituent loadings on 
the land and the rainfall within an area, can vary over 
space and time within the area and time period being 
simulated, and are typically not known with any preci- 
sion. Chapter 13 reviews some of these loading models 
and their limitations. 

4. Water Quality Model Processes 

Water quality models can be applied to many different 
types of water system, including streams, rivers, lakes, 
reservoirs, estuaries, coastal waters and oceans. The models 
describe the main water quality processes, and typically 
require the hydrological and constituent inputs (the water 
flows or volumes and the pollutant loadings). These 
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models include terms for dispersive and/or advective trans- 
port depending on the hydrological and hydrodynamic 
characteristics of the water body, and terms for the biolog- 
ical, chemical and physical reactions among constituents. 
Advective transport dominates in flowing rivers. Dispersion 
is the predominant transport phenomenon in estuaries 
subject to tidal action. Lake-water quality prediction is 
complicated by the influence of random wind directions 
and velocities that often affect surface mixing, currents and 
stratification. For this and other reasons, obtaining reliable 
quality predictions for lakes is often more difficult than for 
streams, rivers and estuaries. In coastal waters and oceans, 
large-scale flow patterns and tide are the most important 
transport mechanisms. 

The development and application of water quality 
models is both a science and an art. Each model reflects 
the creativity of its developer, the particular water quality 
management problems and issues being addressed, the 
available data for model parameter calibration and 
verification, the time available for modelling and associ- 
ated uncertainty, and other considerations. The fact that 
most, if not all, water quality models cannot accurately 
predict what actually happens does not detract from their 
value. Even relatively simple models can help managers 
understand the real world prototype and estimate at least 
the relative, if not actual, change in water quality associ- 
ated with given changes in the inputs resulting from 
management policies or practices. 

4.1. Mass-Balance Principles 

The basic principle of water quality models is that of mass 
balance. A water system can be divided into different 
segments or volume elements, also called ‘computational 
cells’. For each segment or cell, there must be a mass 
balance for each water quality constituent over time. Most 
water quality simulation models simulate quality over a 
consecutive series of discrete time periods, At. Time is 
divided into discrete intervals t and the flows are assumed 
constant within each of those time period intervals. For 
each segment and each time period, the mass balance of 
a substance in a segment can be defined. Components of 
the mass balance for a segment include: first, changes by 
transport (Tr) into and out of the segment; second, changes 
by physical or chemical processes (P) occurring within the 


segment; and third, changes by sources/discharges to or 
from the segment (S). 


M‘ + At = M[ + At 


AMj 

TT 


+ At 


hr 


AM; 

At 


+ At 


"A Mi 

r At 


The mass balance has the following components: 


the mass in computational cell i at the beginning of a 
time step t: M- 

the mass in computational cell i at the end of a time 
step t: M‘ + At 

changes in computational cell i 
r AM, ^ 


by transport: 


At 


Jjr 


changes in computational cell i by physical, 

( AM t ^ 

(bio)chemical or biological processes: 

changes in computational cell i by sources 

f AM^ 


(e.g. wasteloads, river discharges): 


At 


A 


Changes by transport include both advective and disper- 
sive transport. Advective transport is transport by flowing 
water. Dispersive transport results from concentration 
differences. Dispersion in the vertical direction is impor- 
tant if the water column is stratified, and dispersion in the 
horizontal direction can be in one or two dimensions. 
Dispersion, as defined here, differs from the physical 
concept of molecular diffusion as it stands for all 
transport that is not advective. 

Changes by processes include physical processes such 
as re-aeration and settling, (bio)chemical processes 
such as adsorption, transformation and denitrification, 
and biological processes such as primary production and 
predation on phytoplankton. Water quality processes 
convert one substance to another. 

Changes by sources include the addition of mass by 
wasteloads and the extraction of mass by intakes. Mass 
entering over the model boundaries can be considered 
a source as well. The water flowing into or flowing out 
of the modelled segment or volume element (the compu- 
tational cell) is derived from a water quantity (possibly 
hydrodynamic) model. 
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To model the transport of substances over space, a 
water system is divided in small segments or volume 
elements. The complete ensemble of all the segments or 
elements is called the grid or s chematization. Each compu- 
tational cell is defined by its volume and its dimensions in 
one, two or three directions (Ax, Ay, A z) depending on 
the nature of the schematization (ID, 2D or 3D). Note 
that the cell dimensions Ax, Ay and Az do not have to be 
equal. The computational cell can have any rectangular 
shape. A computational cell can share surface areas with 
other cells, the atmosphere, or the bottom sediment or 
coast line. 

The following sections will look at the transport 
processes in more detail, defining parameters or variables 
and their units in terms of mass M, length L and time T. 


4.1.1. Advective Transport 

The advective transport, T A (M/T), of a constituent at a site 
x 0 is the product of the average water velocity, v Xq (L/T), at 
that site, the surface or cross-sectional area, A (L 2 ), through 
which advection takes place at that site, and the average 
concentration, C Xq (M/L 3 ), of the constituent: 

T a = v x XAXC* (12.3) 

X 0 X 0 X 0 V 7 


4.1.2. Dispersive Transport 


The dispersive transport, (M/T), across a surface area 
is assumed to be proportional to the concentration 


gradient 


dC 

dx 


at site x 0 times the surface area A. 

x=x 0 


Letting D Xq (L 2 /T), be the dispersion or diffusion coeffi- 
cient at site x 0 : 



X AX 


dC 

dx 


(12.4) 


Dispersion is done according to Fick’s diffusion law. The 
minus sign originates from the fact that dispersion causes 
net transport from higher to lower concentrations, and so 
in the opposite direction of the concentration gradient. 
The concentration gradient is the difference of concentra- 
tions per unit length, over a very small distance across the 
cross section: 


dC _ Yiyyi ^*+0.5Ax ^x— 0.5Ax 

dx x Ax— >0 Ax 


(12.5) 


Dispersion coefficients should be calibrated or be 
obtained from calculations using turbulence models. 


4.1.3. Mass Transport by Advection and Dispersion 


If the advective and dispersive terms are added and the 
terms at a second surface at site x 0 + Ax are included, 
a one dimensional equation results: 


M\ + At =M\ + A f X 


V x„ Cxx T\' + Ax Cc +Ax 


„ dC 

Dx °dx 


dC 

\ 

+ d Xo+Ax — 

X 0 

x 0 +Ax J 


X A 


( 12 . 6 ) 


or equivalently: 
M[ +At = M[ + At 


X 


Q.x C x Q x +AxCx n + Ax H.X Ax 


dC 

dx 


dC 

' f^x 0 +Ax^-x 0 +Ax ^ 


(12.7) 


x 0 +Ax y 


where Q Xq (L 3 /T) is the flow at site x 0 . 

If the previous equation is divided by the volume and 
the time interval At, then the following equation results in 
one dimension: 


A t s-^t 

C A 


ac 

° X ° +Ax dx 


+Ax 


dC 

Dx °dx 


At 


Ax 


+ 


V xXx 0 


^Xq+Ax^Xq+Ax 


Ax 


( 12 . 8 ) 


Taking the asymptotic limit At — > 0 and Ax — > 0, the 
advection-diffusion equation for one dimension results: 


dC 

dt 


dx 


D- 


d£ 

dx 


— (vC) 

dx 


(12.9) 


The finite volume method for transport is a computa- 
tional method of solving the advection-diffusion 
equation. The accuracy of the method will be related to 
the size of Ax, A (A = AyAz) and At. 
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By adding terms for transport in the y and z-direction, 
a three-dimensional model is obtained. Taking the 
asymptotic limit again will lead to a three-dimensional 
advection-diffusion equation 

dC_ d 2 C 3C d 2 C dC 

x d^~ Vx X y dX~ Vy X 

+ D z —y + V z — + S + f R (C,t ) (12.10) 

dz oz 

with dispersion coefficients Dj defined for each direction. 
If source terms ‘S’ and ‘/ R ’ are added as shown in the equa- 
tion above, the so-called advection-diffusion reaction 
equation emerges. The additional terms represent: 

• Discharges or ‘wasteloads’ (S): these source terms are 
additional inflows of water or mass. As many source 
terms as required may be added to Equation 12.10. 
These could include small rivers, discharges of indus- 
tries, sewage treatment plants, small wasteload outfalls 
and so on. 

• Reaction terms or ‘processes’ (f R ). 

Processes can be split into physical and other processes. 
Examples of physical processes are: 

• settling of suspended particulate matter 

• water movement not affecting substances, like 
evaporation 

• volatilization of the substance itself at the water 
surface. 

Examples of other processes are: 

• biochemical conversions like ammonia and oxygen 
forming nitrite 

• growth of algae (primary production) 

• predation by other animals 

• chemical reactions. 

These processes are described in more detail in the 
remaining parts of this section. 

The expression D(3C/3X)— vC in Equation 12.9, 
multiplied by the area A, is termed the total flux (M/T). 
Flux due to dispersion, DA(3C/dX), is assumed to be 
proportional to the concentration gradient over distance. 
Constituents are transferred by dispersion from higher 
concentration zones to lower ones. The coefficient of 
dispersion D (L 2 /T) depends on the amplitude and 
frequency of the tide, if applicable, as well as upon the 


turbulence of the water body. It is common practice to 
include in this dispersion parameter everything affecting 
the distribution of C other than advection. The term vAC 
is the advective flux caused by the movement of water 
containing the constituent concentration C (M/L 3 ) at a 
velocity rate v (L/T) across a cross-sectional area A (L 2 ). 

The relative importance of dispersion and advection 
depends on the degree of detail with which the velocity 
field is defined. A good spatial and temporal description 
of the velocity field within which the constituent is being 
distributed will reduce the importance of the dispersion 
term. Less precise descriptions of the velocity field, such 
as averaging across irregular cross sections or approxi- 
mating transients by steady flows, may lead to a 
dominance of the dispersion term. 

Many of the reactions affecting the decrease or increase 
of constituent concentrations are often represented 
by first-order kinetics that assume the reaction rates are 
proportional to the constituent concentration. While 
higher-order kinetics may be more correct in certain 
situations, predictions of constituent concentrations 
based on first-order kinetics have often been found to be 
acceptable for natural aquatic systems. 


4.2. Steady-State Models 

A steady state means no change over time. If we consider 
a water body, for example a river, this means there are no 
changes in the concentrations with time. In this case the 
left hand side of Equation 12.9, dC/dt , equals 0. Assume 
the only sink is the natural decay of the constituent 
defined as kC where k, (T -1 ), is the decay rate coefficient 
or constant. Now Equation 12.9 becomes 

0 = D d 2 C/dX 2 - vdC/dX - kC (12.11) 

Equation 12.11 can be integrated, since river reach 
parameters A, D, k, v, and Q are assumed constant. For a 
constant loading, W c (M/T) at site X = 0, the concentra- 
tion C at any distance X will equal 

C(X) = (W c /Qm) exp[(v/2D)(l + m)X] X < 0 

(W c /Qm) exp[(v/2D)(l - m)X] X > 0 (12.12) 

where 

m = (1 + (4feD/v 2 )) 172 (12.13) 
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Note from Equation 12.12 that the parameter m is always 
equal to or greater than 1 and that the exponent of e is 
always negative. Hence, as the distance X increases in 
magnitude, either in the positive or negative direction, the 
concentration C(X) will decrease if k > 0. The maximum 
concentration C occurs at X = 0 and is W c /Qm. 

C( 0) = W c /Qm (12.14) 


These equations are plotted in Figure 12.4. 

In flowing rivers not under the influence of tidal 
actions the dispersion is usually small. Assuming the 
dispersion coefficient D is 0, the parameter m defined by 
Equation 12.13 is 1. Hence, when D = 0, the maximum 
concentration at X = 0 is W c /Q. 

C(0) = W c /Q if D = 0. (12.15) 


Assuming D = 0 and v, Q and k > 0, Equation 12.12 
becomes 


C(X) = 


'0 

XW C /Q ) exp[— kX/v] 


X<0 

X>0 


(12.16) 


The above equation for X > 0 can be derived from 
Equations 12.12 and 12.13 by noting that the term (1 — m) 
equals (1 — m)( 1 + m)/( 1 + m) = (1 — m 2 )/ 2 when D is 0. 
Thus when D is 0 the expression v/2D)(l — m)X in Equation 
12.12 becomes — kXIv. The termX/v is sometimes denoted 
as a single variable representing the time of flow: the time 
flow Q takes to travel from site X = 0 to some other down- 
stream site for a distance of X. 

As rivers approach the sea, the dispersion coefficient D 
increases and the net downstream velocity v decreases. 
Because the flow Q equals the cross-sectional area A times 
the velocity v, Q = Av, and since the parameter m can 
be defined as (v 2 + 4feD) 1/2 /v, then as the velocity 
v approaches 0, the term Qm = Av(v 2 + 4feD) 1/2 /v 
approaches 2 A(feD) 172 . The exponent vX(l ± m)/2D in 
Equation 12.12 approaches ±X(fe/D) 1/2 . 

Hence for small velocities, Equation 12.12 becomes 

_ r (W c /2 A(JeD) 1/2 ) exp[+ X(Je/D) 1/2 ] X< 0 
i ( W c /2 A(kD) 1/2 ) exp [- X(fc/D) 1/2 ] X > 0 

(12.17) 

Here, dispersion is much more important than advective 
transport and the concentration profile approaches a 
symmetric distribution, as shown in Figure 12.4, about 
the point of discharge at X = 0. 



Figure 12.4. Constituent concentration distribution along 
a river or estuary resulting from a constant discharge of that 
constituent at a single point source in that river or estuary. 


Water quality management models are often used to 
assess the effect of pollutant loadings on ambient waters 
and to compare the results with specific water quality 
standards. The above steady-state equations can be 
used to construct such a model for estimating the 
wastewater removal efficiencies required at each waste- 
water discharge site that will result in an ambient 
stream quality that meets the standards along a stream 
or river. 

Figure 12.5 shows a schematic of a river into which 
wastewater containing constituent C is being discharged 
at four sites. Assume that maximum allowable concentra- 
tions of the constituent C are specified at each of those 
discharge sites. To estimate the necessary reduction 
in these discharges, the river must be divided into 
approximately homogenous reaches. Each reach can be 
characterized by constant values of the cross-sectional 
area, A, dispersion coefficient, D, constituent decay rate 
constant, k, and velocity, v, associated with some ‘design’ 
flow and temperature conditions. These parameter values 
and the length, X, of each reach can differ; hence, the 
subscript index i will be used to denote the particular 
parameter values for the particular reach. These reaches 
are shown in Figure 12.5. 
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Figure 12.5. Optimization model 
for finding constituent removal 
efficiencies, Ri, at each 
discharge site / that result in 
meeting stream quality 
standards, C/ max , at least total 
cost. 



minimize XjCost;(7?;) 
subject to: 
mass balances: 

C l = W l (l-R 1 )KQm) l +C 2 TB l 

C 2 = [C\TF i(Qm)i + C 3 TB 2 (Qm) 2 ] / {Qm) 2 

C 3 = fVK 3 ri-/e 3 ) + C 2 TF 2 (Qm) 2 + C 7 TB 3 (Qm) 3 \ / (Qm) 3 

C 4 = [W 4 (l-R 4 )/(Qm) 4 + C 5 TB 4 

C 5 = [C 4 TF 4 (Qm ) 4 + C 6 TB 5 (Qm) 5 ] / (Qm) s 

Cj = [C 5 TF 5 (Qm) 5 + C 7 TB 6 (Qm) 6 ] / ( Qm )j 

C-, = [C 3 TF 3 (Qm) 3+ C 6 TF 6 (Qm) 6 + CBgTB 1 (Qm) 7 ] / (Qm)-j 

r„ = rw 8 n -/? + CnTFJOmX 1 / fftnl. 


• quality standards and maximum removal 
efficiencies: 

C/ < <T ax Vi ; < /q max Vi 



In Figure 12.5 each variable Q represents the 
constituent concentration at the beginning of reach i. 
The flows Q represent the design flow conditions. For 
each reach i the product (Q{m^) is represented by (Qm) t . 
The downstream (forward) transfer coefficient, TF b 
equals the applicable part of Equation 12.12, 

TF i = exp[(v/2D)(l - m)X] (12.18) 

as does the upstream (backward) transfer coefficient, TB t . 

TB i = exp[(v/2D)(l + m)X] (12.19) 

The parameter m is defined by Equation 12.13. 

Solving a model such as the one shown in Figure 12.5 
does not mean that the least-cost wasteload allocation 
plan will be implemented, but least cost solutions can 
identify the additional costs of other imposed constraints, 
for example, to ensure equity or extra safety. Models like 
this can be used to identify the cost-quality tradeoffs 
inherent in any water quality management programme. 
Other than economic objectives can also be used to obtain 
other tradeoffs. 

The model in Figure 12.5 incorporates both advection 
and dispersion. If upstream dispersion under design 
streamflow conditions is not significant in some reaches, 
then the upstream (backward) transfer coefficients, TB b 
for those reaches i will equal 0. 


4.3. Design Streamflows for Water Quality 

In streams and rivers, the water quality may vary signifi- 
cantly, depending on the water flow. If wasteload 
discharges are fairly constant, a high flow serves to dilute 
the waste concentration, while where there is a low 
flow concentrations may become undesirably high. It is 
therefore common practice to pick a low-flow condition 
for judging whether or not ambient water quality 
standards are being met. This can also be seen from 
Equations 12.12, 12.14, 12.15, and 12.16. This often is 
the basis for the assumption that the smaller (or more 
critical) the design flow, the more likely it is that the 
stream quality standards will be met. This is not always 
the case, however. 

Different regions of the world use different design 
low- flow conditions. One example of such a design flow, 
which is used in parts of North America, is the minimum 
seven-day average flow expected once in ten years on 
average. Each year the lowest seven-day average flow is 
determined, as shown in Figure 12.6. The sum of each of 
the 365 sequences of seven average daily flows is divided 
by seven, and the minimum value is selected. This is the 
minimum annual average seven-day flow. 

These minimum seven-day average flows for each 
year of record define a probability distribution whose 
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minimum annual average minimum 

7-day flow 7-day flow period 


J I I 



Figure 12.6. Portion of annual flow time series showing 
low flows and the calculation of average seven and 
fourteen-day flows. 


cumulative probabilities can be plotted. As illustrated 
in Figure 12.7, the particular flow on the cumulative 
distribution that has a 90% chance of being exceeded 
is the design flow. It is the minimum annual average 
seven-day flow expected once in ten years. This flow is 
commonly called the 7Q10 flow. Analyses have shown 
that this daily design flow is exceeded about 99% of the 
time in regions where it is used (NRC, 2001). This means 
that there is on average only a 1% chance that any daily 
flow will be less than this 7Q10 flow. 

Consider now any one of the river reaches shown in 
Figure 12.5. Assume an initial amount of constituent 
mass, M, exists at the beginning of the reach. As the reach 
volume, QAt, increases due to the inflow of less polluted 
water, the initial concentration, M/QAt, will decrease. 
However, the flow velocity will increase, and thus the 
time, At, it takes to transport the constituent mass to the 
end of that reach will decrease. This means less time for 
the decay of the constituent. Thus, wasteload allocations 
that meet ambient water quality standards during low- 
flow conditions may not meet them under higher-flow 
conditions, which are observed much more frequently. 
Figure 12.8 illustrates how this might happen. This does 
not suggest low flows should not be considered when 
allocating wasteloads, but rather that a simulation of 
water quality concentrations over varying flow conditions 
may show that higher-flow conditions at some sites 
are even more critical and more frequent than are the 
low- flow conditions. 

Figure 12.8 shows that for a fixed mass of pollutant at 
X = 0, under low- flow conditions the more restrictive 



Figure 12.7. Determining the minimum seven-day annual 
average flow expected once in ten years, designated 
7Q10, from the cumulative probability distribution of annual 
minimum seven-day average flows. 



Figure 12.8. Increasing the streamflows decreases initial 
concentrations but may increase downstream concentrations. 

(lower) maximum pollutant concentration standard in the 
downstream portion of the river is met, but that same 
standard is violated under more frequent higher-flow 
conditions. 

4.4. Temperature 

Temperature affects almost all water quality processes 
taking place in water bodies. For this reason, it may 
be important to model temperature when it may vary 
substantially over the period of interest, or when the 
discharge of heat into water bodies is to be managed. 
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Temperature models are based on a heat balance in the 
water body. A heat balance takes into account the sources 
and sinks of heat. The main sources of heat in a water 
body are short-wave solar radiation, long-wave atmos- 
pheric radiation, conduction of heat from the atmosphere 
to the water and direct heat inputs. The main sinks of heat 
are long-wave radiation emitted by the water, evaporation 
and conduction from the water to the atmosphere. 
Unfortunately, a model with all the sources and sinks of 
heat requires measurements of a number of variables and 
coefficients that are not always readily available. 

One temperature predictor is the simplified model that 
assumes an equilibrium temperature T e (°C) will be 
reached under steady-state meteorological conditions. 
The temperature mass balance in a volume segment 
depends on the water density p (g/cm 3 ), the heat capacity 
of water, c p (cal/g /°C), and the water depth h (cm). The net 
heat input, K H (T e — T) (cal/cm 2 /day), is assumed to be 
proportional to the difference of the actual temperature, 
T, and the equilibrium temperature, T e (°C). 

dT/dt = K H (T e - T)/pc p h (12.20) 

The overall heat exchange coefficient, K H (cal/cm 2 /day/°C), 
is determined in units of Watts/m 2 /°C (1 cal/cm 2 /day/°C 
= 0.4840 Watts/m 2 /°C) from empirical relationships that 
include wind velocity, dew point temperature and actual 
temperature T (°C) (Thomann and Mueller, 1987). 

The equilibrium temperature, T e , is obtained from 
another empirical relationship involving the overall heat 
exchange coefficient, K H , the dew point temperature, T d , 
and the short-wave solar radiation, H s (cal/cm 2 /day), 

T e = T d + (H s /K h ) (12.21) 

This model simplifies the mathematical relationships of a 
complete heat balance and requires less data. 

4.5. Sources and Sinks 

Sources and sinks of pollutants include wasteloads, and 
the physical and biochemical processes that alter those 
wasteloads. External inputs of each constituent would have 
the form W/QAt or W/(A X AX) where W (M/T) is the loading 
rate of the constituent and QAt or A x AX (L 3 ) represents 
the volume of water into which the mass of waste W is 
discharged. Constituent growth and decay processes are 
discussed in the remaining parts of this Section 4. 


4.6. First-Order Constituents 

The first-order models are commonly used to predict 
water quality constituent decay or growth. They can 
represent constituent reactions such as decay or growth in 
situations where the time rate of change (dC/dt) in the 
concentration C of the constituent, say organic matter 
that creates a biochemical oxygen demand (BOD), is 
proportional to the concentration of either the same or 
another constituent concentration. The temperature- 
dependent proportionality constant k c (1/day) is called 
a rate coefficient or constant. In general, if the rate of 
change in some constituent concentration Cj is propor- 
tional to the concentration C t , of constituent i, then 

dC/dt = a ij ?e j 6>[ T_20) C i (12.22) 

where 6 { is the temperature correction coefficient for k { at 
20 °C and T is the temperature in °C. The parameter a p is 
the grams of Cj produced (a p > 0) or consumed (a tJ - < 0) 
per gram Q. For the prediction of BOD concentration 
over time, C i — Cj = BOD and a p — a B0D = — 1 in 
Equation 12.22. Conservative substances, such as salt, 
will have a decay rate constant k of 0. The concentration 
of conservative substances depends only on the amount 
of water, that is, dilution. 

The typical values for the rate coefficients k c and 
temperature coefficients of some constituents C are 
in Table 12.1. For bacteria, the first-order decay rate (k B ) 
can also be expressed in terms of the time to reach 
90% mortality (t 90 , days). The relationship between these 
coefficients is given by k B = 2.3/t go . 

4.7. Dissolved Oxygen 

Dissolved oxygen (DO) concentration is a common 
indicator of the health of the aquatic ecosystem. DO was 
originally modelled in the Ohio River (US) by Streeter 
and Phelps (1925). Since then a number of modifica- 
tions and extensions of the model have been made 
relating to the number of sinks and sources of DO being 
considered, and how processes involving the nitrogen 
cycle and phytoplankton are being modelled, as 
illustrated in Figure 12.9. 

The sources of DO in a water body include re-aeration 
from the atmosphere, photosynthetic oxygen production 
from aquatic plants, denitrification and DO inputs. The 
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Table 12.1. Typical values of the first-order 
decay rate, k, and the temperature correction 
factor, 6, for some constituents. 


total coliform bacteria (freshwater) 

1. 0-5.5 

- a 

l/day 

total coliform bacteria (sediments) 

0.14-0.21 

- a 

l/day 

total coliform bacteria (seawater) 

0.7-3. 0 

- a 

l/day 

fecal coliform bacteria (seawater) 

37-1 10 

- a 

l/day 

BOD (no treatment) 

0.3-0.4 

- a 

l/day 

BOD (activated sludge treatment) 

0.05-0. 1 

- a 

l/day 

carbofuran 

0.03 

- b 

l/day 

DDT 

0 . 0 - 0.10 

- b 

l/day 

PCB 

0.0-0.007 

- b 

l/day 

pentachlorophenol 

0.0-33.6 

- b 

l/day 


constituent 


rate constant k 


units 


constituent 

0 


units 

coliform bacteria (freshwater) 

1.07 

- b 

— 

coliform bacteria (saltwater) 

1 . 10 

- b 

— 

BOD 

a - Thomann and Mueller ( 1 987) 

1 .04 - a 

b - Schnoor (1996) 




Figure 12.9. The dissolved 
oxygen interactions in a water 
body, showing the decay 
(satisfaction) of carbonaceous, 
nitrogenous and sediment 
oxygen demands and water body 
re-aeration or deaeration (if 
supersaturation occurs at the 
air-water interface). 
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sinks include oxidation of carbonaceous and nitrogenous 
material, sediment oxygen demand and respiration by 
aquatic plants. 

A0 2 /A t = loads + transport + re-aeration 

+ net primary production + denitrification 
— mineralization —nitrification — SOD 

The rate of re-aeration is assumed to be proportional to 
the difference between the saturation concentration, DO sat 
(mg/1), and the concentration of dissolved oxygen, DO 
(mg/1). The proportionality coefficient is the re-aeration 
rate k r (1/day), defined at temperature T = 20 °C , which 
can be corrected for any temperature T with the 
coefficient 6^ 20 \ The value of this temperature correc- 
tion coefficient, 0, depends on the mixing condition of 
the water body. Values generally range from 1.005 to 
1.030. In practice, a value of 1.024 is often used (Thomann 
and Mueller, 1987). The re-aeration rate constant is a 
sensitive parameter. There have been numerous equations 
developed to define this rate constant. Table 12.2 lists 
some of them. 

The saturation concentration, DO sat , of oxygen in 
water is a function of the water temperature and salinity 
(chloride concentration, Cl (g/m 3 )), and can be approxi- 
mated by 

DO sat = {14.652 - 0.41022T + (0.089392 T) 2 

- (0.042685T) 3 }{1 - (Cl/100,000)} (12.23a) 

Elmore and Hayes (1960) derived an analytical expres- 
sion for the DO saturation concentration, DO sat (mg/1), as 
a function of temperature (T, °C): 


DO sat = 14.652 - 0.41022T 4- 0.007991T 2 

- 0.000077774T 3 (12.23b) 

Fitting a second-order polynomial curve to the data 
presented in Chapra (1997) results in: 

DO sat = 14.407 - 0.3369 T + 0.0035 T 2 (12.23c) 

as is shown in Figure 12.10. 

Because photosynthesis occurs during daylight hours, 
photosynthetic oxygen production follows a cyclic, 
diurnal, pattern in water. During the day, oxygen 
concentrations in water are high and can even become 
supersaturated, i.e. concentrations exceeding the satura- 
tion concentration. At night, the concentrations drop due 
to respiration and other oxygen-consuming processes. 



Figure 12.10. Fitted curve to the saturation dissolved oxygen 
concentration (mg/l) as a function of temperature (°C). 


Table 12.2. Some equations for 
defining the re-aeration rate 
constant, k r (day -1 ). 


w water and wind velocity (m/s) water depth (m) 

‘E 

3 

k r = mass transport coefficient for reaeration (m/day) / (water depth) 

= 5.026 (water velocity) 0 - 969 / (water depth) 1 673 (Churchill, 1962) 

= 3.95 (water velocity) 0 - 5 / (water depth) 1 - 5 (O'Connor and Dobbiens, 1 958) 

= (scale factor) 3.95 (water velocity) 0 - 5 / (water depth) 1 - 5 

= 5.344 (water velocity) 0 - 670 / (water depth) 185 (Owens, Edwards, Gibb, 1964) 

= 5.13 (water velocity) / (water depth) 1333 (Langbien, Durum, 1 967) 

= { 0.065 (wind velocity) 2 + 3.86 [(water velocity) / (water depth)] 0 5 } / (water depth) 

(van Pagee 1 978, Delvigne 1 980) 
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One can distinguish between the biochemical oxygen 
demand from carbonaceous organic matter (CBOD, mg/1) 
in the water, and that from nitrogenous organic matter 
(. NBOD , mg/1) in the water. There is also the oxygen 
demand from carbonaceous and nitrogenous organic 
matter in the sediments (SOD, mg/l/day). These oxygen 
demands are typically modelled as first-order decay 
reactions with decay rate constants k CB0D (1/day) for 
CBOD and k NB0D (1/day) for NBOD. These rate constants 
vary with temperature, hence they are typically defined 
for 20 °C. The decay rates are corrected for temperatures 
other than 20 °C using temperature coefficients 6 CB0D 
and 0 NBOD , respectively. 

The sediment oxygen demand SOD (mg/l/day) is usually 
expressed as a zero-order reaction, that is, a constant 
demand. One important feature in modelling NBOD is 
ensuring the appropriate time lag between when it is 
discharged into a water body and when the oxygen demand 
is observed. This lag is in part a function of the level of 
treatment in the wastewater treatment plant. 

The dissolved oxygen (DO) model with CBOD, NBOD 
and SOD is 

dDO/dt = — ^cwh^cbod^CBOD 
~ kmoD dg~ 0 \ fNBOD 

+ fe r 6» r (T_20) (DO sat - DO) - SOD (12.24) 
dCBOD/dt = - k CBOD 0^ T B oD 0) CBOD (12.25) 

dNBOD/df = - fe NBOD 6>^oD 0) NBOD (12.26) 

The mean and range values for coefficients included in 
these dissolved oxygen models are shown in Table 12.3. 

4.8. Nutrients and Eutrophication 

Eutrophication is the progressive process of nutrient 
enrichment of water systems. The increase in nutrients 
leads to an increase in the productivity of the water 
system, which may result in an excessive increase in 
the biomass of algae or other primary producers, such as 
macrophytes or duck weed. When it is visible on the 
surface of the water it is called an algae bloom. Excessive 
algal biomass could affect the water quality, especially 
if it causes anaerobic conditions and thus impairs the 
drinking, recreational and ecological uses. 

The eutrophication component of the model relates 
the concentration of nutrients and the algal biomass. 
Taking the example shown in Figure 12.11, consider 


the growth of algae A (mg/1), depending on phosphate 
phosphorus, P (mg/1), and nitrite/nitrate nitrogen, N n mg/1), 
as the limiting nutrients. There could be other limiting 
nutrients or other conditions as well, but here consider only 
these two. If either of these two nutrients is absent, the 
algae cannot grow, regardless of the abundance of the other 
nutrient. The uptake of the more abundant nutrient will 
not occur. 

To account for this, algal growth is commonly 
modelled as a Michaelis-Menten multiplicative effect; in 
other words, the nutrients have a synergistic effect. Model 
parameters include a maximum algal growth rate jl 
(1/day) times the fraction of a day, f d , that rate applies 
(Figure 12.12), the half saturation constants K P and K N 
(mg/1) (shown as K c in Figure 12.13) for phosphate and 
nitrate, respectively, and a combined algal respiration and 
specific death rate constant e (1/day) that creates an 
oxygen demand. The uptake of phosphate, ammonia and 
nitrite/nitrate by algae is assumed to occur in proportion 
to their contents in the algae biomass. Define these 
proportions as a P , a A and a N , respectively. 

In addition to the above parameters, one needs to 
know the amounts of oxygen consumed in the oxidation 
of organic phosphorus, P 0 , and the amounts of oxygen 
produced by photosynthesis and consumed by respira- 
tion. In the model below, some average values have 
been assumed. Also assumed are constant temperature 
correction factors for all processes pertaining to any 
individual constituent. This reduces the number of 
parameters needed, but is not necessarily realistic. Clearly 
other processes as well as other parameters could be 
added, but the purpose here is to illustrate how these 
models are developed. Users of water quality simulation 
programs will appreciate the many different assumptions 
that can be made and the large amount of parameters 
associated with most of them. 

The source and sink terms of the relatively simple 
eutrophication model shown in Figure 12.11 can be 
written as follows: 


For algae biomass: 

dA/dt = jufA T - 20) [P/(P + Kp)}[N n /(N n + K n )] 


X A - e6>i r_20) A 

(12.27) 

For organic phosphorus: 


dP /d t = - k 0 (T ~ 20) P 

va± 0 /vai, op w op 1 0 

(12.28) 
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Table 12.3. Typical values of 
parameters used in the 
dissolved oxygen models. 


parameter 


value 


units 


kr, slow, deep rivers 0. 1 -0.4 

- a 

l/day 

kr, typical conditions 0.4- 1.5 

- a 

l/day 

kr, swift, deep rivers 1 .5-4.0 

- a 

l/day 

k r , swift, shallow rivers 4.0- 1 0.0 

- a 

l/day 

kcBOD , untreated discharges 0.35 (0.20-0.50) 

-b 

l/day 

kcBOD , primary treatment 0.20 (0.10-0.30) 

-b 

l/day 

kcBOD, activated sludge 0.075 (0.05-0.10) 

-b 

l/day 

0C8OD 1-04 

- a 

— 

1.047 

- a 

— 

1.04 (1.02-1.09) 

- c 

— 

0 r 1.024 (1.005-1.030) 

- c 

___ 


sediment oxygen demand * 


value 


units 


municipal sludge (outfall vicinity) 

4 

(2-10) 

- c d 

g / m 2 / day 

municipal sewage sludge 

1.5 

(1-2) 

- c d 

g 0 2 / m 2 / day 

sandy bottom 

0.5 

(0.2- 1.0) 

- c d 

g 0 2 / m 2 / day 

mineral soils 

0.07 

(0.05-0.1) 

- c d 

g 0 2 / m 2 / day 

natural to low pollution 

0.1-10.0 

- a 

g 0 2 / m 2 / day 

moderate to heavy pollution 

5-10 


- a 

g 0 2 / m 2 / day 


a - Schnoor (1 996) c - Thomann and Mueller (1 987) 

b - Chapra ( 1 997) d - Bowie et al. ( 1 985) 

* value has to be divided by the water height (m) 


For phosphate phosphorus: 
dP/dt = - juj d ei T ~ 20) [P/(P + Kp)] 
X[N n /(N„ + K N )la P A 
For organic nitrogen: 
dN /d t = - k 0 (r-2O) N 

on w on Ly o 

For ammonia-nitrogen: 

dN a /d£ = - Jufd0i T ~ 2O) mP + Kp)] 

X [N n /(N„ + K N )]a A A 
+ k 0 (t ~ 2O) N - k 0 (t ~ 2O) N 

1 '^on^on iN o K aYa Ly a 


For nitrate-nitrogen: 


(12.29) 

(12.30) 


(12.31) 


dN n /d£ = - juf d 0A~ 2O) [P/(P + K p )] 

X [N„/(N n + K N )] flN A 
+ k 0 (t ~ 2O) N - k 0 (t ~ 2O) N 

1 ' t fl u fl Ly a '^n w n ly n 


For dissolved oxygen: 


dDO/dt — k CBOD 9^ 0 ^CBOD 


- 4.57kX T ~ 20) K - 2 k op 0%- 20) P o 
+ (1.5 fif d - 2 e)6>f- 20) A 
+ k r 0 r (T ~ 2O) (DO 5al - DO) - SOD 


(12.32) 


(12.33) 
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Figure 12.11. The dissolved oxygen, 
nitrogen and phosphorus cycles, and 
phytoplankton interactions in a water 
body, showing the decay (satisfaction) of 
carbonaceous and sediment oxygen 
demands, re-aeration or deaeration of 
oxygen at the air-water interface, 
ammonification of organic nitrogen in the 
detritus, nitrification (oxidation) of 
ammonium to nitrate-nitrogen and 
oxidation of organic phosphorus in the 
sediment or bottom layer to phosphate 
phosphorus, phytoplankton production 
from nitrate and phosphate consumption, 
and phytoplankton respiration and death 
contributing to the organic nitrogen and 
phosphorus. 



Figure 12.12. Calculation of the fraction, f d , of the maximum 
growth rate constant, ju, to use in the algal growth equations. 
The fraction f d is the ratio of actual production zone/potential 
production zone. 

Representative values of the coefficients for this model are 
shown in Table 12.4. 

Because the growth of phytoplankton cannot occur 
without nutrients, the eutrophication modelling must be 
coupled with that of nutrients. Nutrient modelling must 
include all the different biochemical forms of the nutri- 
ents, primarily nitrogen and phosphorus, as well as all the 
interactions between the different forms. The sum of all 
these interactions is referred to a ‘nutrient cycling’. 



Figure 12.13. Defining the half saturation constant for a 
Michaelis-Menten model of algae. The actual growth rate 
constant = juC/[C + K c ). 


The nitrogen cycle includes ammonium (NH 4 -N) and 
nitrate/nitrite (represented as N0 3 -N) as the main forms 
of dissolved nitrogen in water. Furthermore, nitrogen is 
present in algae, as well as detritus, resulting from algae 
mortality, and suspended (non-detritus) organic nitrogen. 
Nitrogen can also be present in different forms in 
the bottom sediment. Two important reactions in the 
nitrogen nutrient cycle are nitrification and denitrifica- 
tion, which affect the flux of ammonium and nitrate in the 
water column. 
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Table 12.4. Typical values of 
coefficients in the eutrophication 
model. 


parameter 


value 


units 


k N 

half saturation 1 0-20 

- a 

figN/l 


50-200 

- c 

MS N0 3 /l 


10 (1-20) 

-b 

MgN/l 

kp 

half saturation 1 -5 

- a 

Mg P/I 


20-70 

- c 

Mg P/I 


10 

-b 

Mg P/I 

Op 

stoichoimetric ratio 0.0 1 2-0.0 1 5 

- c 

mg P/mg A 

a 

N 

stoichoimetric ratio 0.08-0.09 

- c 

mg N0 3 /mg A 

M 

maximum algae growth rate 1.5 (1 .0-2.0) 

-b 

l/day 

e 

death algae rate 0.2-8 

- c 

l/day 


0. 1 (0.05-0.025) 

-b 

l/day 


a - Thomann and Mueller ( 1 987) 
b - Schnoor (1996) 
c - Bowie et al. (1985) 


Nitrification is the conversion of ammonium to nitrite 
and finally nitrate, requiring the presence of oxygen: 

NH 4 + 20 2 -> NO 3 + 2H 2 0 + 2H + (12.34) 

Denitrification is the process occurring during the break- 
down (oxidation) of organic matter by which nitrate is 
transformed to nitrogen gas, which is then usually lost 
from the water system. Denitrification occurs in anaerobic 
conditions: 

NO 3 -> N 2 (12.35) 

The phosphorus cycle is simpler than the nitrogen 
cycle because there are fewer forms in which phosphorus 
can be present. There is only one form of dissolved phos- 
phorus, orthophosphorus (also called orthophosphate, 
PO4-P). Like nitrogen, phosphorus also exists in algae, 
in detritus and other organic material, as well as in the 
bottom sediment. Unlike nitrogen, there can also be 
inorganic phosphorus in the particulate phase. 

Further details of the nutrient cycles are given in 
Section 5. 

4.9. Toxic Chemicals 

Toxic chemicals, also referred to as ‘micro-pollutants’, 
are substances that at low concentrations can impair the 
reproduction and growth of organisms, including fish 


and human beings. These substances include heavy 
metals, many synthetic organic compounds (organic 
micro-pollutants) and radioactive substances. 

4.9.1. Adsorbed and Dissolved Pollutants 

An important characteristic of many of these substances 
is their affinity with the surface areas of suspended or 
bottom sediments. Many chemicals preferentially sorb 
onto surfaces of particulate matter rather than remaining 
dissolved in water. To model the transport and fate of 
these substances, the adsorption-desorption process, 
estimations of the suspended sediment concentration, 
resuspension from the bottom and settling are required. 

Figure 12.14 depicts the adsorption-desorption and 
first-order decay processes for toxic chemicals and their 
interaction in water and sediment. This applies to the 
water and sediment phases in both the water body and 
the bottom sediments. 

The adsorption-desorption model assumes (conve- 
niently but not always precisely) that an equilibrium 
exists between the dissolved (in water) and adsorbed 
(on sediments) concentrations of a toxic constituent 
such as a heavy metal or organic contaminant. This 
equilibrium follows a linear relationship. The slope of that 
linear relation is the partition coefficient K p (litres/kg). 
This is shown in Figure 12.15. 
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waterbody 




Figure 12.14. Schematic of the adsorption/desorption and 
decay processes of various toxic chemicals in water bodies 
and bottom sediments. 



Figure 12.15. Defining the partition coefficient K p (litres/kg) 
as the slope of the fixed ratio between concentrations of a 
constituent in the water and sediment phases of either a 
water body or bottom sediments. Different constituents have 
different partition coefficients, when they apply. 


Each partition coefficient K p (mg/kgDW/mg/1 water or 
1/kg) is defined as the ratio of the particulate concentra- 
tion Cp of a micro-pollutant (mg/kgDW or mg/kgC) 
divided by the dissolved concentration Q' of a micro- 
pollutant (mg/1 water). 

K p = C'/Q (12.36) 

Representative values of partition coefficients K p are given 
in Table 12.5. 

The presence of a micro-pollutant in a water system is 
described by the total concentration (sum of dissolved 
and particulate concentrations), the total particulate 
concentration and the total dissolved concentration for 
each water and sediment compartment. The particulate 
and dissolved concentrations are derived from the total 
concentration and the respective fractions. 

Because the fate of most micro-pollutants is largely 
determined by adsorption to particulate matter, sus- 
pended inorganic and organic matter (including phyto- 
plankton) has to be included in the model in most cases. 
It may be necessary to include dissolved organic matter 
as well. 

The adsorbed fractions in the water column are subject 
to settling. The fractions in the sediment are subject to 
resuspension. The adsorbed fractions in the sediment 
can also be removed from the modelled part of the water 
system by burial. 

The rates of settling and resuspension of micro- 
pollutants are proportional to the rates for particulate 
matter. An additional process called bioturbation leads to 


K p parameter 


value 


units 


arsenic 

I0 4 

l/kg 

heavy metals (Cd,Cu,Cr,Zn) 

I0 4 -I0 6 

l/kg 

benzo(a)pyrene 

I0 4 -I0 5 

l/kg 

lead 

I0 5 -I0 6 

l/kg 

PCB 

I0 5 -I0 6 

l/kg 

plutonium-239 

o 

o 

l/kg 

methoxychlor 

I0 4 

l/kg 

napthalene 

I0 3 

l/kg 


Table 12.5. Typical values of partition coefficients in toxic 
chemical model from Thomann and Mueller (1987). 


the redistribution of the micro-pollutant among sediment 
layers. Bioturbation is caused by the physical activity of 
organisms, and affects both the particulate and dissolved 
phases but at different rates. Bioturbation is taken into 
account by means of dispersion coefficients. 

For modelling purposes, it is important to know how 
much of a toxic chemical is present as a dissolved 
constituent, as opposed to adsorbed. Assuming partition 
coefficients apply to a particular toxic constituent, the 
concentration, C w , of that constituent in the water body is 
divided into a dissolved fraction (f dw ) and an adsorbed 
fraction (f aw ); 

C w = (f dw +f a JC w (12.37) 

In turn, the adsorbed fraction is composed the fractions 
of a micro-pollutant adsorbed to inorganic particulate 
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matter, /im, dead particulate organic matter, fpoc, and 
algae, falg. The total micro-pollutant concentration, 
C w (mg/m 3 ) is the sum of all these fractions: 

C w = (fdw + Jim 4- fpoc 4- falg)C w (12.38) 

Considering the simple division into dissolved and 
adsorbed fractions (f dw and / aw ), these fractions depend 
on the partition coefficient, K p , and on the suspended 
sediment concentration, SS (mg/1). The proportions of 
the total constituent concentration in the water body, C w , 
dissolved in the water, DC W (mg/1), and adsorbed to the 
suspended sediments, AC W (mg/1) are defined as 

DC w = f dw C w (12.39) 

AC w =f aw C w (12.40) 

where the fractions 

f dw = 1/(1 + K p SS) (12.41) 

fa w = 1 - Cw = K p SS/(l + KpSS ) (12.42) 

Similarly, in the bottom sediments, the dissolved concen- 
tration DC S (mg/1) and adsorbed concentration AC S (mg/1) 
are fractions, and / fls , of the total concentration C s (mg/1); 

DC s =f ds C s (12.43) 

AC s =f as C s (12.44) 

These fractions are dependent on the sediment porosity, 
0, and density, p s (kg/1): 

fds = l/[0 + Ad - 0) K p] (12.45) 

fas= 1 -/* = ([0 + AU - 0Kp] - 1)/ 

[</) + A(1 - 0Kp] (12.46) 

First order decay occurs in the water and sediment phases 
only in the dissolved fraction with decay rate constants k w 
and k s (1/day), respectively. Thus: 

dC w /dt = - fe w C“ 20) CwQ 

- f aw C wS + fasQr (12.47) 

dC s /dt = - k s e s (T - 20) f ds c s 

+ faw C wS - fas C s r (12.48) 

In the above two equations the parameter s represents 

the mass of settling sediments (mg/day), r the mass of 
resuspension sediments (mg/day), and 6 the temperature 
correction coefficient of the constituent at temperature 


T = 20 °C. If data are not available to distinguish 
between the values of the decay rate constants k in water 
and on sediments, they may be assumed to be the same. 
Similarly, for the values of the temperature correction 
coefficients 6. The settling and resuspension of suspended 
solids can be determined each day from a sediment 
model. 

4.9.2. Heavy Metals 

The behaviour of heavy metals in the environment 
depends on their inherent chemical properties. Heavy 
metals can be divided into different categories accord- 
ing to their dissolved form and redox status. Some 
metals, including copper, cadmium, lead, mercury, 
nickel, tin and zinc, form free or complexed cations 
when dissolved in water, e.g. Cu 2+ or CuCl - . The 
soluble complexes are formed with negatively charged 
ions such as chlorine, oxygen or dissolved organic 
compounds. These heavy metals also tend to form 
poorly soluble sulphides under chemically reducing 
conditions. These sulphides generally settle in bottom 
sediments and are essentially unavailable ecologically. 
Other metals such as arsenic and vanadium are present 
as anions in dissolved form. The differences between 
groups of metals have important consequences for the 
partitioning of the metals among several dissolved and 
particulate phases. 

Metals are non-decaying substances. The fate of heavy 
metals in a water system is determined primarily by 
partitioning to water and particulate matter (including 
phytoplankton), and by transport. The partitioning 
divides the total amount of a pollutant into a ‘dissolved’ 
fraction and several ‘adsorbed’ fractions (as described in 
Equations 12.39-12.42). The fractions of a metal that are 
adsorbed onto particulate matter are influenced by all the 
processes that affect particulate matter, such as settling 
and resuspension. 

Partitioning is described in general by sorption to 
particulates, precipitation in minerals, and complexation 
in solution. Complexation with inorganic and organic 
ligands can be considered explicitly in connection with 
the other processes. Sorption can be modelled as an 
equilibrium process (equilibrium partitioning) or as the 
resultant of slow adsorption and desorption reactions 
(kinetic formulations). In the latter case, partitioning 
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is assumed to proceed at a finite rate proportional 
to the difference between the actual state and the 
equilibrium state. 

To describe the fate of certain heavy metals in reducing 
environments, such as sediment layers, the formation of 
metal sulphides or hydroxides can be modelled. The 
soluble metal concentration is determined on the basis of 
the relevant solubility product. The excess metal is stored 
in a precipitated metal fraction. 

Sorption and precipitation affect the dissolved metal 
concentration in different ways. Both the adsorbed 
and dissolved fractions increase at increasing total 
concentration as long as no solubility product is 
exceeded. When it is, precipitation occurs. 

4.9.3. Organic Micro-pollutants 

Organic micro-pollutants are generally biocides (such as 
pesticides or herbicides), solvents or combustion 
products, and include substances such as hexachloro- 
hexane, hexachlorobenzene, PCBs or polychloro- 
biphenyls, benzo-a-pyrene and fluoranthene (PAHs or 
polycyclic aromatic hydrocarbons), diuron and linuron, 
atrazine and simazine, mevinfos and dichlorvos, and 
dinoseb. 

The short-term fate of organic micro-pollutants in a 
water system is determined primarily by partitioning 
to water and organic particulate matter (including 
phytoplankton), and by transport. Additional processes 
such as volatilization and degradation influence organic 
micro-pollutant concentrations (this is in contrast to 
heavy metals, which do not decay). Many toxic organic 
compounds have decay (or ‘daughter’) products that 
are equally, if not more, toxic than the original 
compound. The rates of these processes are concentration 
and temperature dependent. 

Organic micro-pollutants are generally very poorly 
soluble in water and prefer to adsorb to particulate 
matter in the water, especially particulate organic matter 
and algae. Therefore, the fractions of a micro-pollutant 
adsorbed to inorganic matter, /im, dead particulate 
organic matter, jpoc, the dissolved fraction of a micro- 
pollutant, fd, and algae, falg , add up to the total 
micro-pollutant concentration, C (mg/m 3 ); 

C = (fd + fim + jpoc + falg)C (12.49) 


The fractions are functions of the partition coefficients 
K p (for algae (m 3 /gC), for inorganic matter (m 3 /gDW) 
and for dead particulate organic matter (m 3 /gC)), the 
individual concentrations C (for algae biomass (gC/m 3 ), for 
inorganic matter (gDW/m 3 ) and for dead particulate organic 
matter (gC/m 3 )), and the porosity 0(m 3 water/m 3 bulk). In 
surface water the value for porosity is 1 . 

fd = </>/[</) + Kpalg Calg + Kpim Cim 

+ Kppoc Cpoc] (12.50) 

fim = (1 — fd) Kpim Cim/ 

[Kpalg Calg + Kpim Cim + Kppoc Cpoc] (12.51) 

fpoc = (1 — fd) Kppoc Cpoc/ 

[Kpalg Calg + Kpim Cim + Kppoc Cpoc] (12.52) 

falg = (1 — fd — fim — fpoc) (12.53) 

In terms of bulk measures, each partition coefficient K p 
(see Equation 12.36) also equals the porosity (j) times the 
bulk particulate concentration Cp (mg/m 3 bulk) divided 
by the product of the dissolved (mg/1 bulk) and particu- 
late (mg/m 3 bulk) bulk concentrations, Cd Cs, all times 
10 6 mg/kg. 

K p = 10 6 (j)Cp/(CdCs) (12.54) 

Partitioning can be simulated based on the above equilib- 
rium approach or according to slow sorption kinetics. For 
the latter, the rate, dCp/dt, of adsorption or desorption 
(mg/m 3 /day) depends on a first order kinetic constant 
ksorp (day -1 ) for adsorption and desorption times 
the difference between equilibrium particulate concentra- 
tion Cpe of a micro-pollutant (mg/m 3 bulk) and the 
actual particulate concentration Cp (mg/m 3 bulk) of a 
micro-pollutant. 

dCp/dt = ksorp(Cpe — Cp) (12.55) 

The kinetic constant for sorption is not temperature 
dependent. All other kinetic constants for micro-pollutants 
are temperature dependent. 

Mass-balance equations are similar for all micro- 
pollutants except for the loss processes. 

Metals are conservative substances that can be 
transformed into various species either through complexa- 
tion, adsorption or precipitation. Organic micro-pollutants 
are lost by volatilization, biodegradation, photolysis, 
hydrolysis and overall degradation. Most of these processes 
are usually modelled as first-order processes, with associ- 
ated rate constants. 
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Volatilization is formulated according to the double 
film theory. The volatization rate dCd/dt (mg/m 3 /day) 
of dissolved micro-pollutant concentrations, Cd (mg/m 3 
water), in water depends on an overall transfer coefficient, 
kvol (m/day), for volatilization and the depth of the water 
column, H (m); 

dCd/dt = - kvol Cd/H (12.56) 

The numerator ( kvol Cd) is the volatilization mass flux 
(mg/m 2 /day). 

This equation is only valid when the atmospheric 
concentration is negligibly small, which is the normal 
situation. 

All other loss rates such as biodegradation, photolysis, 
hydrolysis or overall degradation (mg/m 3 /day) are usually 
modelled as 

dC/dt = ~kC (12.57) 

where C is the total concentration of a micro-pollutant 
(mg/m 3 ), and k is a (pseudo) temperature-dependent 
first-order kinetic rate constant for biodegradation, 
photolysis, hydrolysis or overall degradation (day -1 ). 

4.9.4. Radioactive Substances 

The fate of most radionuclides, such as isotopes of iodine 
( 131 1) and cesium ( 137 Cs), in water is determined prima- 
rily by partitioning to water and particulate matter 
(including phytoplankton), by transport and by decay. 
Cesium (Ce + ) adsorbs to particulate inorganic matter, to 
dead particulate organic material, and to phytoplankton, 
both reversibly and irreversibly. The irreversible fraction 
increases over time as the reversible fraction gradually 
transforms into the irreversible fraction. Radioactive 
decay proceeds equally for all fractions. Precipitation of 
cesium does not occur at low concentrations in natural 
water systems. 

Iodine is only present in soluble form as an anion 
(10 J) and does not adsorb to particulate matter. 
Consequently, with respect to transport, iodine is only 
subject to advection and dispersion. 

Concentrations of radionuclides, C R (mg/m 3 ), are 
essentially conservative in a chemical sense, but they 
decay by falling apart in other nuclides and various types 
of radiation. The radioactive decay rate (mg/m 3 /day) is 
usually modelled as a first order process involving a 
kinetic radioactive decay constant, kdec (day -1 ). This 


kinetic constant is derived from the half-life time of the 
radionuclide. The initial concentration may be expressed 
as radioactivity, in order to simulate the activity instead of 
the concentration. These state variables can be converted 
into each other using: 

Ac = 10 -3 N A kdecC R /[ 86400 Mw\ (12.58) 

where 

Ac — activity of the radionuclide (Bq/m 3 /s) 

N a = Avogadro’s number (6.02 * 10 23 mole) 

Mw = molecular weight of the radionuclide (g/mole). 

4.10. Sediments 

As discussed in the previous section, sediments play 
an important role in the transport and fate of chemical 
pollutants in water. Natural waters can contain a mixture 
of particles ranging from gravel (2 mm to 20 mm) or 
sand (0.07 mm to 2 mm) down to very small particles 
classified as silt or clay (smaller than 0.07 mm). The very 
fine fractions can be carried as colloidal suspension, for 
which electrochemical forces play a predominant role. 
Considering its large adsorbing capacities, the fine 
fraction is characterized as cohesive sediment. Cohesive 
sediment can include silt and clay particles as well as 
particulate organic matter such as detritus and other 
forms of organic carbon, diatoms and other algae. Since 
flocculation and adsorbing capacities are of minor 
importance for larger particles, they are classified as 
non-cohesive sediment. 

The behaviour of this fine-grained suspended matter 
affects water quality in several ways. First, turbidity and 
its effect on the underwater light climate is an important 
environmental condition for algae growth. The presence 
of suspended sediment increases the attenuation of 
light in the water column, which leads to an inhibition 
of photosynthetic activity, and hence a reduction in 
primary production. Second, the fate of contaminants in 
waters is closely related to suspended solids due to 
their large adsorbing capacities. Like dissolved matter, 
sediment is transported by advection and by turbulent 
motion. In addition, the fate of the suspended cohesive 
sediment is determined by settling and deposition, as well 
as by bed processes of consolidation, bioturbation and 
resuspension. 
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4.10.1. Processes and Modelling Assumptions 

From a water quality perspective, the three most important 
sediment processes are sedimentation, resuspension or 
‘erosion’, and ‘burial’ of sediment. 

The modelling of sedimentation and erosion processes 
originates in part from the Partheniades-Krone concept 
(Partheniades, 1962; Krone, 1962). In this concept, the 
bottom shear stress plays an essential role in defining 
whether or not sedimentation of suspended particles or 
erosion of bed material will occur. Sedimentation takes 
place when the bottom shear stress drops below a critical 
value. Resuspension occurs when the bottom shear stress 
exceeds a critical value. The calculation of bed shear stress 
is discussed in Section 4.10.5. 

For purposes of modelling chemical concentrations 
in the bottom sediment, one or more sediment layers 
are usually defined. Each sediment layer is assumed 
homogeneous (well mixed). The density of the layer can 
vary, depending on the variable sediment layer composi- 
tion. The porosity within a given layer is assumed constant 
and user-defined. Water quality models do not generally 
consider horizontal transport of bed sediments. This hori- 
zontal transport would result in a change in the amount of 
sediment present in the bottom, and hence a change in the 
thickness of the layer of non-cohesive sediment. 

4.10.2. Sedimentation 

A characteristic feature of cohesive sediments is their 
ability to form aggregates of floes that settle to the bottom, 
in a process called sedimentation. Whether a particle will 
settle to the bottom depends on its size and density, and 
the chemical conditions of the surrounding water system. 
Various laboratory and field measurements show that the 
suspended matter concentration strongly influences the 
aggregation process and thereby the settling velocities of 
the aggregates (Krone, 1962). Strong floes are denser and 
have larger settling velocities. 

Sediment floe aggregation depends on the chemical 
and physical properties of the sediment, and on salinity 
and turbulence. At high sediment concentrations (several 
g/1) the particles hinder each other, decreasing their 
settling velocity. Turbulence affects the flocculation and 
therefore the settling velocity in two opposing ways. 
Increased turbulence will increase the collisions between 


particles, resulting in larger floes with higher settling 
velocities, but increased turbulent shear stresses will 
break up the floes and decrease their settling velocities. 

Sedimentation occurs when the bottom shear stress 
drops below a critical value. Total shear stress is the sum of 
the shear stresses induced by flows and wind waves. 
Sedimentation rates can depend on salinity concentrations. 

For any sediment particle j that settles out of the water 
column, the rate of settling depends on the flow velocity 
shear stress, T (kg/m/sec 2 ), at the bottom surface-water 
interface, the critical shear stress, Tj, for the substance j, 
a zero order sedimentation rate, ZSj (g/m 2 /day), and the 
sedimentation velocity, VSj (m/day), for the settling 
substance, Cj (g dry wt./m 3 /day), all divided by the depth, 
H (m), of settling in the water column. 

If the shear stress Tat the bottom- water interface is less 
than the critical shear stress Tj for the substance j, settling 
is assumed to take place. The rate of decrease of the 
substance (g dry wt./m 3 /day) in the water column due to 
settling is 

dCj/dt = - [1 - (T/Tj)] C ZSj + VSfp/H T< Tj 

(12.59) 

The corresponding flux (g dry wt./m 2 of sediment layer 
area/day) of substance j onto the sediment layer (accre- 
tion) is 

dCf ed /dt = - (dCj/dt) H (12.60) 

Sedimentation always results in an increase of the 
substance settling in the upper sediment bed layer. The 
quantity of sedimentation in one model time-step cannot 
exceed the available amount of substance in the water 
column. 

No net resuspension is assumed to take place if 
sedimentation occurs. 

4.10.3. Resuspension 

If the shear stress T exceeds the critical shear stress Tj for 
the type substance j, then no net settling takes place. 
Instead, resuspension can take place (also called 
‘erosion’). 

Erosion of cohesive bed material occurs when the bed 
shear forces exceed the resistance of the bed sediment. The 
resistance of the bed is characterized by a certain critical 
erosive strength (bottom shear stress). This critical stress is 
determined by several factors, such as the chemical 
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composition of the bed material, particle size distribution 
and bioturbation. Erosion of sediment is induced by the 
bed stress due to flow velocities, tidal and wind-induced 
advective flows, and surface waves. Erosion is directly 
proportional to the excess of the applied shear stress over 
the critical erosive bottom shear stress. One formula for 
erosion of homogeneous beds is based on Partheniades 
(1962). The erosion/resuspension flux is limited by the 
available amount of sediment on the sea bed. Typically a 
one-layer homogeneous bed is assumed. 

The rate of resuspension (g dry wt./m 2 of bottom layer 
area/day) of substance j going into the water column 
depends on these flow velocity shear stresses, a zero-order 
resuspension flux ZR (g/m 2 /day) if any, a resuspension 
rate constant kRj (day -1 ), and the amount of substance 
C- ed in the top sediment layer (g/m 2 ); 

dCj/dt = [( T/Tj) - 1 KZRj + kRjCj Sed )/H Tj 

(12.61) 

The corresponding flux (g dry wt./m 2 sediment layer 
area/day) of substance Cj sed from the sediment layer area 
(erosion) is 

dq ed /dt = - (dCj/dt) H (12.62) 

No sedimentation is assumed if resuspension occurs. 
Resuspension results in a decrease of sediment in the 
upper sediment bed layer. The sedimentation in one 
model time step cannot exceed the available amount 
of substance in the sediment layer. The values of ZRj and 
T strongly depend on the sediment properties and envi- 
ronmental parameter values. 

For water quality considerations, the chemical compo- 
sition of the resuspending sediment is assumed to be 
the same as that of the bottom sediment. Resuspension of 
chemical component concentrations, Nj (g/m 3 /day), in 
a substance from the sediment layer is simply the 
substance, C sed resuspension flux (g dry wt./m 2 /day) 
(Equation 12.62) times the fraction,/^, of component j in 
the substance divided by the depth, H, of water column: 

dN/dt = frj(dC sed /dt)/H (12.63) 


4.10.4. Burial 

‘Burial is a term used to convey that bed sediment is 
no longer available for resuspension because it has 
been covered by much newer sediment, i.e. it has been 


buried. This concept is important especially with respect 
to contaminated sediments, as it implies that buried 
contaminated sediments are isolated from the overlying 
water column and no longer pose a threat to the water 
quality. 

Burial is a sink for sediments that are otherwise 
susceptible for resuspension. Assuming a fixed active 
sediment layer thickness, burial occurs when this thick- 
ness is exceeded. Establishing this requires knowledge 
of the depth, hj (m/g dry wt/m 2 ) of particle size 
substance j per unit dry wt per square metre. Hence, if 
the total depth of the active sediment layer, XjCf ed hj, is 
greater than the assumed maximum depth, D max (m), of 
active sediment layer, then burial will result. The 
decrease in available sediment substance j will be in 
proportion to its contribution, Cf ed hj , to the total depth, 
Xj Cj Sed hj, of the sediment layer. The burial flux will 
equal a burial rate constant k B (day -1 and usually 
set equal to 1) times the amount Cf ed of substance 
j times the depth ratio times the excess depth ratio, 



(12.64) 


4.10.5. Bed Shear Stress 

The bed shear stress directly influences the sedimentation 
and erosion rates. It depends on the flow (currents) and 
the wind-generated (and sometimes human- generated) 
surface waves. For sedimentation-erosion processes, 
it is usually assumed that the total bed shear stress, T 
(Newton/m 2 or kg/m/s 2 ) due to flow, 7" fLow , and waves, 
T wave , are additive: 

^ ~ Tlow Tvave (12.65) 

The bed shear stress T flow for depth-averaged flow 
depends on the water density (1,000 kg/m 3 ), the hori- 
zontal flow velocity, U h (m/s), acceleration of gravity, 
g (9.81 m/s 2 ), and a Chezy coefficient, Cz (m°' 5 /s). 

T now = 1,000 (9.81 )Ul/Cz 2 (12.66) 
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The Chezy coefficient, Cz (m° 5 /s), is defined as a function 
of water depth, H (m), and bottom roughness length 
(Nikuradse equivalent roughness length) (m), Rough (m), 
or Manning’s roughness coefficient, n (s/m 173 ): 

Cz = 18 log 10 (12 H/Rough ) (12.67) 

Cz = (H) 1/6 /n (12.68) 

Both Nikuradse’s roughness coefficient and Manning’s 
roughness coefficient may be changing due to bed load 
movements. Here (and in most models) they are assumed 
fixed. 

Surface waves are caused by wind stress on the water 
surface. The magnitude of the waves depends on the 
wind conditions, wind duration, water depth and bottom 
friction. Wave fields are commonly described by the 
significant wave height, significant wave period and 
wavelength. Waves induce a vertical circular movement 
(orbital velocity) that decreases with depth. The waves 
exert friction forces on the bed during propagation. 

The magnitude of the bed shear stress, T wave , due to 
waves depends on a wave friction factor, / w , the density of 
water, (1,000 kg/m 3 ) and the effective orbital horizontal 
velocity at the bed surface, U 0 . 

T wave = 0.25 (1,000) / w l/ 0 2 (12.69) 

The wave friction factor, / w , and the effective orbital 
horizontal velocity at the bed surface, U 0 , are functions 
of three wave parameters: the significant wave height, 
H s (m), the mean wave period, T m (s), and mean wave 
length, L m (m) (see also Figure 12.16). Also required is 
the depth of water, H (m). 

The effective horizontal bottom velocity due to waves 


is defined as 

U 0 = 7rH s /[T m sinh(2;rH/L m )] (12.70) 

The friction or shear factor, / w , can be calculated in two 
ways (Monbaliu, et al., 1999). One way is 

f w = 0.16 [ Rough/(U 0 T m /2x)l 05 (12.71) 

The other way uses a factor depending on a parameter A 
defined as 

A = H s /[2 Rough sinh(2;rH/L m )] (12.72) 

If A > 1.47 then 

f w = exp{— 5.977 + 5.123H -0194 } (12.73) 

Otherwise 

f w = 0.32 (12.74) 



Figure 12.16. Wave dimensions of significant wave height, 
period and length. 


4.11. Lakes and Reservoirs 

The water quality modelling principles discussed above 
are applicable to all different types of water systems: 
streams, rivers, lakes, estuaries and even coastal or 
ocean waters. This section presents some of the unique 
aspects of water quality modelling in lakes. The physi- 
cal character and water quality of rivers draining 
into lakes and reservoirs are governed in part by the 
velocity and volume of river water. The characteristics 
of the river water typically undergo significant changes 
as the water enters the lake or reservoir, primarily 
because its velocity reduces. Sediment and other 
material that were carried in the faster-flowing water 
settle out in the basin. 

The structure of the biological communities also 
changes from organisms suited to living in flowing 
waters to those that thrive in standing or pooled waters. 
There are greater opportunities for the growth of 
algae (phytoplankton) and the development of 
eutrophication. 

Reservoirs typically receive larger inputs of water (as 
well as soil and other materials carried in rivers) than 
lakes do, and as a result, may receive larger pollutant 
loads. However, because of greater water inflows, flushing 
rates are more rapid than in lakes. Thus, although reser- 
voirs may receive greater pollutant loads than lakes, they 
have the potential to flush out the pollutants more 
rapidly. Reservoirs may therefore exhibit fewer or less 
severe negative water quality or biological impacts than 
lakes for the same pollutant load. 

The water quality of lakes and reservoirs is defined by 
variables measured within the water basin. Although 
there are many variables of limnological significance, 
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water quality is typically characterized on the basis of 
conditions such as: 

• water clarity or transparency (greater water clarity 
usually indicates better water quality) 

• concentration of nutrients (lower concentrations 
indicate better water quality) 

• quantity of algae (lower levels indicate better water 
quality) 

• oxygen concentration (higher concentrations are 
preferred for fisheries) 

• concentration of dissolved minerals (lower values 
indicate better water quality) 

• acidity (a neutral pH of 7 is preferred). 

Many waste chemical compounds from industry, some 
with toxic or deleterious effects on humans and/or other 
water-dependent organisms and products, are discharged 
into lakes and reservoirs. Some can kill aquatic organisms 
and damage irrigated crops. Inadequate water purification 
resulting in the discharge of bacteria, viruses and other 


organisms into natural waters can be a primary cause of 
waterborne disease. Although dangerous to human health 
worldwide, such problems are particularly severe in 
developing countries. 

There can be major differences between deep and 
shallow lakes or reservoirs. Deep lakes, particularly in 
non-tropical regions, usually have poorer water quality in 
their lower layers, due to stratification (see Section 4.11 .3). 
Shallow lakes do not exhibit this depth differentiation 
in quality. Their more shallow, shoreline areas have 
relatively poorer water quality because those sites are 
where pollutant inputs are discharged and have a greater 
potential for disturbance of bottom muds and the like. 
The water quality of a natural lake usually improves as 
one moves from the shoreline to the deeper central part. 

In contrast, the deepest end of a reservoir is usually 
immediately upstream of the dam. Water quality usually 
improves along the length of a reservoir, from the shallow 
inflow end to the deeper, ‘lake-like’ end near the dam, as 
shown in Figure 12.17. 


Figure 12.17. Longitudinal 
zonation of water quality and 
other variables in reservoirs 
(UNDP, 2000). 



riverine zone 


o narrow, channelized basin 

o relatively high flow rates 

o high suspended solids 
i= low light availability at 
depth 

o nutrient supply by 
advection 

a relatively high nutrient 
levels 

o light-limited primary 
productivity 

o cell losses primarily 
by sedimentation 

o organic matter sypply 
primarily allochthonous 

o more eutrophic 


transitional zone 


o broader, deeper basin 

o reduced flow rates 

o reduced suspended solids 
■=1 light availability at 
depth 

o advective nutrient supply 
reduced 


o primary productivity 
relatively high 

o cell losses by grazing 
and sedimentation 

o intermediate 
o intermediate 


lacustrine zone 


o broad, deep, lake-like basin 

o low flow rates 

o relatively clear 

= high light availability at 
depth 

o nutrient supply by 
internal recycling 
= relatively low nutrient 
levels 

o nutrient-limited primary 
productivity 

o cell losses primarily 
by grazing 

o organic matter sypply 
primarily autochthonous 

o more oligotrophic 
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Reservoirs, particularly the deeper ones, are also 
distinguished from lakes by the presence of a longitudinal 
gradient in physical, chemical and biological water 
quality characteristics from the upstream river end to the 
downstream dam end. Because of this, reservoirs have 
been characterized as comprising three major zones: an 
upstream riverine zone, a downstream lake-like zone at 
the dam end, and a transitional zone separating these two 
(Figure 12.17). The relative size and volume of the three 
zones can vary greatly in a given reservoir. 

4.11.1. Downstream Characteristics 

Constructing a dam can produce dramatic changes in the 
downstream river channel below it. These are quite unlike 
downstream changes from lakes. Because reservoirs act 
as sediment and nutrient traps, the water at the dam end 
of a reservoir is typically of higher quality than water 
entering the reservoir. This higher-quality water subse- 
quently flows into the downstream river channel below 
the dam. This phenomenon is sometimes a problem, 
in that the smaller the quantity of sediments and other 
materials transported in the discharged water, the greater 
the quantity that it can now pick up and transport. 
Because it contains less sediment, the discharged water 
can scour and erode the streambed and banks, picking up 
new sediment as it continues downstream. This scouring 
effect can have significant negative impacts on the flora, 
fauna and biological community structure in the down- 
stream river channel. The removal of sediments from a 
river by reservoirs also has important biological effects, 
particularly on floodplains. 

Many reservoirs, especially those used for drinking 
supplies, have water release or discharge structures located 
at different vertical levels in their dams (Figure 12.18). 
This allows for the withdrawal or discharge of water 
from different layers within the reservoir, known 
as ‘selective withdrawal’. Depending on the quality of 
the water discharged, selective withdrawal can signifi- 
cantly affect water quality within the reservoir itself, as 
well as the chemical composition and temperature of 
the downstream river. Being able to regulate both 
quantities and qualities of the downstream hydrological 
regimes makes it possible to affect both flora and 
fauna, and possibly even the geomorphology of the 
stream or river. 



Figure 12.18. A multiple-outlet reservoir can be better used 
to regulate the temperature and water quality downstream. 


Constructing a reservoir may have significant social 
and economic implications, including the potential for 
stimulating urban and agricultural development adjacent 
to, and below, the reservoir. These activities can have 
both positive and negative impacts on downstream 
water quality, depending on the nature and size of 
development. 

Agriculture is often the leading source of pollution 
in lakes. Healthy lake ecosystems contain nutrients in 
small quantities from natural sources, but extra inputs of 
nutrients (primarily nitrogen and phosphorus) unbalance 
lake ecosystems. When temperature and light conditions 
are favourable, excessive nutrients stimulate population 
explosions of undesirable algae and aquatic weeds. After 
the algae die, they sink to the lake bottom, where bacte- 
ria consume the available dissolved oxygen as they 
decompose the algae. Fish may die and foul odours may 
result if dissolved oxygen is depleted. 

Heavy metals are another major cause of lake quality 
impairment. Our knowledge of this is mainly due to 
the widespread detection of heavy metals in fish tissue 
samples. Since it is difficult to measure mercury in 
ambient water, and since it bioaccumulates in fish tissue, 
fish samples are commonly used to indicate the level of 
contamination. Common sources of mercury involve 
atmospheric transport from power-generating facilities 
and other ‘smoke-stack’ industries. 

In addition to nutrient and metal siltation, enrich- 
ment by organic wastes that deplete oxygen and by 
noxious aquatic plants affect lakes and reservoirs. Often, 
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several pollutants and processes affect a single lake. For 
example, a process such as removal of shoreline vegeta- 
tion may accelerate erosion of sediment and nutrients 
into a lake. Extreme acidity (low pH) resulting from acid 
rain can eliminate fish in isolated lakes. Urban runoff 
and storm sewers, municipal sewage treatment plants, 
and hydrological modifications are also sources of lake 
pollutants. 

The prediction of water quality in surface water 
impoundments is based on mass-balance relationships 
similar to those used to predict water quality concentra- 
tions in streams and estuaries. There are also significant 
problems in predicting the water quality of lakes or reser- 
voirs compared to those of river and estuarine systems. 
One is the increased importance of wind-induced mixing 
processes and thermal stratification. Another, for reser- 
voirs, is the impact of various reservoir-operating policies. 


4.11.2. Lake Quality Models 

Perhaps the simplest way to begin modelling lakes is 
to consider shallow well-mixed constant-volume lakes 
subject to a constant pollutant loading. The flux of any 
constituent concentration, C, in the lake equals the mass 
input of the constituent less the mass output less losses 
due to decay or sedimentation, if any, all divided by the 
lake volume V (m 3 ). Given a constant constituent input 
rate W c (g/day) of a constituent having a net decay and 
sedimentation rate constant k c (day -1 ) into a lake having 
a volume V (m 3 ) and inflow and outflow rate of 
Q (m 3 /day), then the rate of change in the concentration 
C (g/m 3 /day) is 

dC/dt = (1/V) (W c - QC ~ k c CV ) (12.75) 

Integrating this equation yields a predictive expression of 
the concentration C(t) of the constituent at the end of any 
time period t based in part on what the concentration, 
C(t— 1), was at the end of the previous time period, t— 1. 
For a period duration of At days, 

at) = [W C /(Q + k c V)\ 

[1 - exp{— At((Q/V) + fe c )}] 

+ C(t— 1) exp{— A t((Q/V) + k c )} (12.76) 

The equilibrium concentration, Q, can be obtained by 
assuming each concentration is equal in Equation 12.76 
or by setting the rate in Equation 12.75 to 0, or by setting 


C(t — 1) equal to 0 and letting At go to infinity in Equation 
12.76. The net result is 

C e = W C /(Q + k c V) (12.77) 

The time, t a , since the introduction of a mass input W c 
that is required to reach a given fraction a of the equilib- 
rium concentration (i.e. , C(t)/Q = a) is 


t a = ~ V[ln(l - a)]/(Q + K C V) (12.78) 

Similar equations can be developed to estimate the 
concentrations and times associated with a decrease in 
a pollutant concentration. For the perfectly mixed lake 
having an initial constituent concentration C(0), say after 
an accidental spill, and no further additions, the change 
in concentration with respect to time is 

dC/dt = - C (Q + k c V)/V (12.79) 


Integrating this equation, the concentration C(t) is 
at) = C(0) exp{— t((Q/V) + K c ) (12.80) 


In this case, one can solve for the time t a required for the 
constituent to reach a fraction (1 — a) of the initial con- 
centration C(0) (that is, C(t)/C(0) = (1 — a)). The result is 
Equation 12.78. 

Equation 12.76 can be used to form an optimization 
model for determining the wasteload inputs to this well- 
mixed lake that meet water quality standards. Assuming 
that the total of all natural wasteloads W c (0, inflows and 
outflows Q(t), and the maximum allowable constituent 
concentrations in the lake, C(0 max , may vary among 
different within-year periods £, the minimum fraction, X, 
of total waste removal required can be found by solving 
the following linear optimization model: 

Minimize X (12.81) 

Subject to the following mass-balance and constituent 

concentration constraints for each period t: 

at) = [W C (0(1-X)/(Q(0 + k c V)] 

[1 - exp {—A £((Q(£)/V) + k c )}] 

+ at~ 1) exp{— A £((Q(£)/V) + k c )} (12.82) 

at) < c(0 max (12.83) 


If each period t is a within-year period, and if the waste- 
loadings and flows in each year are the same, then no 
initial concentrations need be assumed and a steady state 
solution can be found. This solution will indicate, for the 
loadings W c (t), the fraction X of waste removal that meet 
the quality standards, C(0 max , throughout the year. 
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4.11.3. Stratified Impoundments 

Many deep impoundments become stratified during 
particular times of the year. Vertical temperature 
gradients arise that imply vertical density gradients. The 
depth-dependent density gradients in stratified lakes 
effectively prevent complete vertical mixing. Particularly 
in the summer season, lakes may exhibit two zones: an 
upper volume of warm water called the epilimnion and 
a lower colder volume called the hypolimnion. The transi- 
tion zone or boundary between the two zones is called the 
thermocline (Figure 12.18). 

Because of lake stratification, many models divide the 
depth of water into layers, each of which is assumed to be 
fully mixed. To illustrate this approach without getting 
into too much detail, consider a simple two-layer lake in 
the summer that becomes a one-layer lake in the winter. 
This is illustrated in Figure 12.19 

Discharges of a mass W c of constituent C in a flow 
Q m (t) into the lake in period t have concentrations of 
W c /Q m (0. The concentration in the outflows from the 
summer epilimnion is Q(0 for each period t in the sum- 
mer season. The concentration of the outflows from the 


winter lake as a whole is C(t) for each period t in the 
winter season. The summer time rates of change in the epil- 
imnion constituent concentrations C e (t) and hypolimnion 
concentrations Q(0 depend on the mass inflow, W c (t), and 
outflow, Q(0Q out (0, the net vertical transfer across thermo- 
cline, (v/D T )[C h (t)V h (t) — C e (t) V e (t)] , the settling on sediment 
interface, sH h (t)C h (t ), and the decay, kC e (t): 

d.C e (t)/dt = (l/V e m(W c (t) - C e (0Q out (0) 

+ (v/D T )[C h (t)V h (t) - C e (0V e m - kC e (t) 

(12.84) 

dC h (£)/d£ = - kC h (t) - (v/D t )[Q(0 - C e (t)V e (t)/V h (t)l 

- sH h (t)C h (t) (12.85) 

In the above two equations, V e and V h (m 3 ) are the time- 
dependent volumes of the epilimnion and hypolimnion 
respectively; k (day -1 ) is the temperature-corrected decay 
rate constant; v (m/day) is the net vertical exchange 
velocity that includes effects of vertical dispersion, 
erosion of hypolimnion, and other processes that transfer 
materials across the thermocline of thickness D r (m); s is 
the settling rate velocity (m/day); and H h (t ) is the average 
depth of the hypolimnion (m). 



Figure 12.19. Lake stratification during 
summer and complete mixing during 
winter season. 
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In the winter season the lake is assumed to be fully 
mixed. Thus for all periods t in the winter season the 


initial concentration of a constituent is: 

ao = c e w e (o + c h m h mv e + v h ] ( 12 . 86 ) 

dC(0/dt = (l/V(0){(W c (t) -at) Q out (0) 

- kC e (t) - sH(t)C(t) (12.87) 

At the beginning of the summer season, each epilimnion 
and hypolimnion concentration will be the same: 

C e (t) = at) (12.88) 

C h (t) = at) (12.89) 


5. An Algal Biomass Prediction 
Model 

An alternative approach to modelling the nutrient, 
oxygen and algae parts of an ecological model has 
been implemented in a simulation model developed by 
WL | Delft Hydraulics called DELWAQ-BLOOM (Los, 
1991; Los et al., 1992; Molen et al., 1994; Smits 2001; 
WL | Delft Hydraulics, 2003). This model is used to 
predict algae growth and mortality, oxygen concentra- 
tions and nutrient dynamics. 

The ecological model DELWAQ-BLOOM calculates 
the advection and dispersion of constituents (state 
variables) in the water column, and the water quality and 
ecological processes affecting the concentrations of the 
constituents. It is based on a three-dimensional version of 
the governing Equation 12.10. The focus here will be on 
the source and sink terms in that equation that define the 
water quality and ecological processes, mostly related to 
algae growth and mortality, mineralization of organic 
matter, nutrient uptake and release, and oxygen produc- 
tion and consumption. 

Lor this discussion, consider three nutrient cycles 
(nitrogen, N, phosphorus, P, and silica, Si) and five 
different groups of algae, (phytoplankton (diatoms, flag- 
ellates and dinoflagellates) or macroalgae (‘attached’ 
or ‘suspended’ Ulva )), suspended and bottom detritus, 
oxygen and inorganic phosphorus particulate matter 
in the bottom sediments. Different predefined sets 
of species are available for marine and freshwater 
ecosystems. 


The model processes relating these substances are all 
interrelated. However, for clarity, the processes can 
be grouped into nutrient cycling, algae modelling, and 
oxygen-related processes. 

5.1. Nutrient Cycling 

The DELWAQ-BLOOM model assumes that algae 
consume ammonia and nitrate in the water column. It 
includes the uptake of inorganic nutrients by bottom 
algae, algae-mortality-producing detritus and inorganic 
nutrients, mineralization of detritus in the water col- 
umn producing inorganic nutrients, and mineralization 
of detritus in the bottom that also produces inorganic 
nutrients. Optionally the model includes algae species 
with the ability to take up atomic nitrogen (N-fixation) 
and detritus (mixotrophic growth). The model accounts 
for the settling of suspended detritus and algae and 
inorganic adsorbed phosphorus, resuspension of 
bottom detritus, release of inorganic bottom nutrients 
to the water, burial of bottom detritus, nitrification 
or denitrification depending on the dissolved oxygen 
concentration, and adsorption/desorption of ortho- 
phosphate. 

5.2. Mineralization of Detritus 

The oxidation or mineralization of the nutrients in detri- 
tus ( DetN , DetP, DetSi) and also of detritus carbon ( DetC ) 
reduces detritus concentrations. The mineralization 
process consumes oxygen and produces inorganic nutri- 
ents (NH 4 , P0 4 , and Si0 2 ). The fluxes, dC/dt, for these 
four constituents C (mg/1 or g/m 3 ) are assumed to be 
governed by first order processes whose temperature 
corrected rate constants are k c O c ^ T ~ 20 ^ (1/day). Thus: 

dC/dt = k c e^ T ~ 20) C (12.90) 

This equation applies in the water column as well as 
in the bottom sediments, but the mineralization rate 
constants, k c ^-2°) 

may differ. This rate constant also 
depends on the stochiometric composition of detritus, 
relative to the requirements of the bacteria. The concen- 
trations of these detritus constituents in the bottom are 
usually expressed in grams per square metre of surface 
area divided by the depth of the sediment layer. 
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5.3. Settling of Detritus and Inorganic 
Particulate Phosphorus 

The rate of settling of nutrients in detritus and inorganic 
particulate phosphorus out of the water column and on to 
the bottom is assumed to be proportional to their water 
column concentrations, C. Settling decreases the concen- 
trations of these constituents in the water column. 

dC/dt = - SR c (C)/H (12.91) 

The parameter SR C is the settling velocity (m/day) of 
constituent concentration C and H is the depth (m) of the 
water column. 

5.4. Resuspension of Detritus and Inorganic 
Particulate Phosphorus 

The rates at which nutrients in detritus and inorganic 
particulate phosphorus are resuspended depend on the 
flow velocities and resulting shear stresses at the bed 
surface-water column interface. Below a critical shear 
stress no resuspension occurs. Resuspension increases the 
masses of these constituents in the water column without 
changing its volume, and hence increases their concentra- 
tions in it. For C B representing the concentration (grams 
of dry weight per cubic metre) of resuspended material in 


the active bottom sediment layer, the flux of constituent 
concentration in the water column is 

dC/dt = RR c C b /H (12.92) 

where RR C (m/day) is the velocity of resuspension 
(depending on the flow velocity) and H is the depth of the 
water column. 

5.5. The Nitrogen Cycle 

The nitrogen cycle considers the water-column compo- 
nents of ammonia (NH 4 -N), nitrite and nitrate (represented 
together as N0 3 -N), algae (AIgN), suspended detritus 
(DefN), and suspended (non-detritus) organic nitrogen 
(OON). In the bottom, sediment bottom detritus ( BDetN ) 
and bottom diatoms ( BDiatN ) are considered. Figure 12.20 
shows this nitrogen cycle. 

5.5.1. Nitrification and Denitrification 

Two important reactions in the nitrogen nutrient cycle are 
nitrification and denitrification. These reactions affect the 
flux of ammonia and nitrate in the water column. Given 
sufficient dissolved oxygen and temperature, nitrifying 
bacteria in the water column transform ammonium to 



Figure 12.20. The nitrogen cycle 
processes. 
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nitrite and then nitrate. This can be considered as one 
reaction, 

NH 4 +20 2 -> NOJ + H 2 0 + 2H + 

that occurs at a rate (gN/m 3 /day) of 

NH 4 /di = - k NHl C 4 20 ) NH 4 (12.93) 

Again, is the temperature corrected rate 

constant (1/day), and NH 4 is the concentration of 
nitrogen in ammonium (gN/m 3 ). 

Bacterial activities decrease as temperatures decrease. 
Bacterial activities also require oxygen. The nitrification 
process stops if the dissolved oxygen level drops below 
about 2 mg/1 or if the temperature T is less than approxi- 
mately 5 °C. 

For each gram of nitrogen in ammonium-nitrogen 
NH 4 -N reduced by nitrification, a gram of nitrate-nitrogen 
NO 3 -N is produced, consuming 2 moles (64 grams) 
oxygen per mole (14 grams) of nitrogen (64/14 = 4.57 
grams of oxygen per gram of nitrogen). Nitrification 
occurs only in the water column. 

In surface waters with a low dissolved-oxygen content, 
nitrate can be transformed to free nitrogen by bacterial 
activity as part of the process of mineralizing organic 
material. This denitrification process can be written as: 

‘organic matter’ + 2N0 3 — > N 2 4- C0 2 + H 2 0 (12.94) 

Nitrate is (directly) removed from the system by means of 
denitrification. The reaction proceeds at a rate: 

dN0 3 /dt = — feNO3^NO 3 20) ^ ( ^ ) 3 (12.95) 

where N0 3 is the concentration of nitrate nitrogen 
(gN/m 3 ). 

This process can occur both in the water column and 
the sediment, but in both cases results in a loss of nitrate 
from the water column. Algae also take up nitrate- 
nitrogen. As with nitrification, denitrification decreases 
with temperature. The reaction is assumed to stop below 
about 5 °C. 

5.5.2. Inorganic Nitrogen 

Ammonia is produced when algae die and the cells release 
their contents into the surrounding water (in a process 
called autolysis) and by the mineralization of organic 
nitrogen in the water and bottom sediment. Ammonia is 
converted to nitrate by nitrification. Algae use ammonia 


and nitrate for growth. Different algae prefer either NH 4 or 
N0 3 nitrogen. Upon death they release part of their nitro- 
gen contents as ammonia. The remaining nitrogen of dying 
algae becomes suspended detritus and suspended ‘other 
organic nitrogen’ (OON). The latter degrades at a much 
slower rate. Algae can also settle to the bottom. Some 
macroalgae ( Ulva ) can be fixed to the bottom, unless wind 
and water velocities are high enough to dislodge them. 

Once in the bottom sediment, planktonic algae die and 
release all their nitrogen contents, as ammonium into the 
water column and to organic nitrogen in the sediment. In 
contrast macroalgae and bottom diatoms are attached 
and are subjected to the normal processes of growth, 
mortality and respiration. 

Suspended detritus and organic nitrogen are formed 
upon the death of algae. Detritus is also produced by 
excretion of phyto- and zooplankton and from resuspen- 
sion of organic matter on and in the sediment. The 
detritus concentration in the water column decreases 
by bacterial decay, sedimentation and filtration by 
zooplankton and benthic suspension feeders. 

Bottom detritus is subject to the processes of miner- 
alization, resuspension and burial. Mineralization of 
bottom detritus is assumed to be slower than that of 
suspended detritus. The ammonia produced from 
mineralization is assumed to go directly to ammonia 
in the water phase. Sedimentation from the water 
column and mortality of algae in the bottom increase the 
bottom pool of bottom detritus. The mineralization 
rate depends on the composition of the detritus, i.e. is a 
function of the nitrogen/carbon and phosphorus/carbon 
ratios. 

Nitrogen is removed from the system by means of 
denitrification, a process that occurs under anoxic condi- 
tions. Burial is a process that puts the material in a ‘deep’ 
sediment layer, and effectively removes it from the active 
system. This is the only removal process for the other 
nutrients (P and Si). 

5.6. Phosphorus Cycle 

The phosphorus cycle (Figure 12.21) is a simplified 
version of the nitrogen cycle. There is only one dissolved 
pool: orthophosphorus, and only one removal process: 
burial. However, unlike nitrogen and silica, there is also 
inorganic phosphorus in the particulate phase (AAP). 
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Figure 12.21. The processes 
involved in the phosphorus 
cycle. 


The phosphorus cycle in the water column includes 
orthophosphate (P0 4 ), algae (AlgP), suspended detritus 
(. DetP), suspended (non-detritus) organic phosphorus 
(OOP), inorganic adsorbed (available) P (AAP), and 
inorganic adsorbed (unavailable) P (LAP). In the bottom 
sediment the cycle includes the bottom detritus ( BDetP ) 
and the bottom inorganic adsorbed P ( BAAP ). 

A reaction specific to the phosphorus cycle is 
the adsorbtion/desorption of particulate inorganic phos- 
phorus. Inorganic phosphorus can be present in the 
aquatic environment in a dissolved form and adsorbed to 
inorganic particles, such as calcium or iron. The transition 
from one form into another is not a first order kinetic 
process, yet in many models desorption of inorganic 
phosphorus is assumed to be such. 

5.7. Silica Cycle 

The silica cycle is similar to the phosphorus cycle except 
that there is no adsorption of silica to inorganic 
suspended solids. Silica is only used by diatoms, so 
uptake by algae depends on the presence of diatoms. The 
silica cycle is shown in Figure 12.22. 

The silica cycle in the water column includes dissolved 
silica (Si), diatoms ( Diat ), suspended detritus ( DetSi ), and 
suspended (non-detritus) organic silica ( OOSi ). In the 


bottom sediment the cycle includes the bottom detritus 
(. BDetSi ). 

5.8. Summary of Nutrient Cycles 

The nutrient cycles just described are based on the 
assumption that nutrients can be recycled an infinite 
number of times without any losses other than due to 
transport, chemical adsorption, denitrification and burial. 
This is an over-simplification of the organic part of the 



Figure 12.22. The processes involved in the silica cycle. 
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nutrient cycles. The elementary composition of living 
algae cells is a complicated function of their characteris- 
tics as well as the environmental conditions. Upon dying, 
the algae cell contents are released into the surrounding 
water. As the cell breaks apart, many of the nutrients are 
in a form that makes them instantly available for algae cell 
growth (autolysis). The remaining material consists of 
more or less degradable substances. Most of this material 
is mineralized either in the water or at the bottom, but a 
small portion degrades very slowly if at all. Most of this 
material settles and is ultimately buried. Resuspension 
delays but does not stop this process by which nutrients 
are permanently removed from the water system. 

For simplicity, all possible removal processes are 
lumped into a single term, which is modelled as burial. 
For example, if a nominal value of 0.0025 per day is used, 
this means that 0.25% of the bottom amount is buried 
each day. 

The same formulation is used for all three nutrients. 
Whether or not this is correct depends on the actual 
removal process. If deactivation is mainly burial into 
deeper layers of the sediment, there is no reason to 
distinguish between different nutrients. Other processes 
such as chemical binding, however, may deactivate 
phosphorus, but not nitrogen or silica. 


5.9. Algae Modelling 

Algal processes include primary production, mortality 
(producing detritus and inorganic nutrients) and grazing, 
settling to become bottom algae and the resuspension of 
bottom algae, the mortality of bottom algae to bottom 
detritus, and the burial of bottom algae. 

The basic behaviour of algae in surface water can be 
illustrated by the two diagrams in Figures 12.23 and 
12.24. These show the nutrient and carbon fluxes for 
diatoms and other algae. Diatoms are distinguished from 
other algae in that they need silicate to grow. 


5.9.1. Algae Species Concentrations 

The module BLOOM II computes phytoplankton within 
DELWAQ-BLOOM. It is based upon the principle of 
competition between different species. The basic variables 
of this module are called ‘types’. A type represents the 



Figure 12.23. Modelling of diatoms. 



Figure 12.24. Modelling of other algae besides diatoms. 


physiological state of a species. Usually a distinction 
is made between three different types: an N-type for 
nitrogen limitation, a P- type for phosphorus limitation 
and an E-type for light energy limitation, whose nutrient 
content is high (luxury uptake). The solution algorithm of 
the model considers all potentially limiting factors and 
first selects the one that is most likely to become limiting. 
It then selects the best-adapted type for that limitation. 
The suitability of a type (its fitness) is determined by the 
ratio of its requirement for that particular resource and 
its growth rate in equal proportion. This means that a 
type can become dominant either because it needs a 
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comparatively small amount of a limiting resource (it is 
efficient) or because it grows rapidly (it is opportunistic). 

Now the algorithm considers the next potentially 
limiting factor and again selects the best-adapted phyto- 
plankton type. This procedure is repeated until it is 
impossible to select a new pair of a type and limiting 
factor without violating (that is, over-exhausting) some 
limiting factor. Thus the model seeks the optimum 
solution consisting of n types and n limiting factors. To 
that purpose, it uses linear programming. Typically, 
BLOOM II considers between three and ten representative 
algae species and six to fifteen types. 

As a further refinement BLOOM II takes the existing 
biomasses of all phytoplankton types into account. These 
are the result of growth, loss processes and transport 
during the previous time steps. Thus, the optimization 
algorithm does not start from scratch. 

As they represent different physiological stages of the 
same species, the transition of one type to another is a rapid 
process with a characteristic time step in the order of a day. 
The equations of the model are formulated in a way that 
allows such a rapid shift between types of the same species. 
A transition between different species is a much slower 
process as it depends on mortality and net growth rates. 

It is interesting that the principle just described, by 
which each phytoplankton type maximizes its own 
benefit, effectively means that the total net production of 
the phytoplankton community is maximized. 

Denote each distinct species subtype (from now 
on called type) by the index k. The BLOOM II module 
identifies the optimum concentration of biomass, B h of 
each algae type k that can be supported in the aquatic 
environment characterized by light conditions and nutri- 
ent concentrations. Using the net growth rate constant Pn h 
it will 

Maximize ^Pn k B k (12.96) 

k 

The sum of the biomass concentrations over all types k is 
the total algae biomass concentration. 

For each algae type, the requirements for nitrogen, 
phosphorus and silica (used only by diatoms) are speci- 
fied by coefficients n ik , the fraction of nutrient i per unit 
biomass concentration of algae type k. 

The total readily available concentration, Q (g/m 3 ) of 
each nutrient in the water column equals the amount in 


the total living biomass of algae, Xk(n ik B k ), plus the 
amount incorporated in dead algae, d b plus that dissolved 
in the water, vq. These mass-balance constraints apply for 
each nutrient i. 

X(% B k ) + d i + w, = C, (12.97) 

k 

The unknown concentration variables B k , d { and vq are 
non-negative. All nutrient concentrations Q are the 
modelled (optionally measured) total concentrations and 
are assumed to remain constant throughout the time 
period (typically a day). 

5.9.2. Nutrient Recycling 

A certain amount of each algae type k dies in each time 
step, and this takes nutrients out of the live phytoplank- 
ton pool. Most of it remains in the detritus and other 
organic nutrient pools; a smaller fraction (in the order of 
30%) is directly available to grow new algae because the 
dead cells break apart (autolysis) and are dissolved in the 
water column. Detritus may be removed to the bottom or 
to the dissolved nutrient pools at rates in proportion to its 
concentration. Needed to model this are the mortality 
rate, M k (day -1 ), of algae type k; the fraction, f p ; of dead 
phytoplankton cells that is not immediately released 
when a cell dies; the remineralization rate constant, m i 
(day -1 ), of dead phytoplankton cells; the fraction, n ik , of 
nutrient i per unit biomass concentration of algae type k; 
and the settling rate constant, s (day -1 ), of dead phyto- 
plankton cells. 

The rate of change in the nutrient concentration of 
the dead phytoplankton cells, dd/dt, in the water column 
equals the increase due to mortality less that which 
remineralizes and that which settles to the bottom; 

dd/dt = ^(f p M k n ik B k ) — m — sd t (12.98) 

k 

Both mortality and mineralization rate constants are 
temperature dependent. 

5.9.3. Energy Limitation 

Algae absorb light for photosynthesis. Energy becomes 
limiting through self-shading when the total light absorp- 
tion of algae, called light extinction, exceeds the 
maximum at which primary production is just balanced 
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by respiration and mortality. For each algae type k there 
exists a specific extinction value K™ ax (1/m) at which this 
is the case. The light intensity can also be too high, which 
means the total extinction is too low (photo-inhibition) 
for growth. This specific extinction value is K™. 
The ranges between K™ m and Kjf ax differ for different 
algal types k. Letting K k (m 3 /m/gdry) represent the 
specific light absorbing extinction constant for living 
material of algae type k, the total extinction due to all 
living algae is 

KL = Z(K k B k ) (12.99) 

k 

Added to this must be the extinction caused by dead cells, 
KD, and the background extinction of the water, KW 
(1/m). As defined here, KW includes the extinction by 
inorganic suspended matter and humic substances. 

Kf n <KL + KD + KW< Kjf ax (12.100) 

In an alternate formulation used in DELWAQ-BLOOM, 
the total light extinction coefficient in the water column 
(in m -1 ) is the sum of background extinction (K b , water 
and non-modelled substances) and extinction due to 
detritus (K det ), inorganic material (K im ), dissolved organic 
matter or yellow substances’ (K dom ) and phytoplankton 

(Kphytk 

Extinction = K b + K det 4- K im + K dom + K phyt (12.101) 

The extinction from each component is calculated by 
its concentration multiplied by its specific extinction 
coefficient. 

The ratio between the various extinction components 
is highly site specific. In general it also varies seasonally. 
The amount of dead cells not yet mineralized is, from 
Equation 12.98, Xfe(/pM fe B fe ). Assuming some fraction e d 
(usually between 0.2 and 0.4) of the extinction rate of 
live cells, 

KD = e d XK k fpM k B k (12.102) 

k 

If the total extinction is not within the range for an algae 
type k, then its concentration B k will be zero. To ensure 
that B k is 0 if the total extinction is outside of its extinc- 
tion range, a 0, 1 binary (integer) unknown variable Z fe 
is needed for each algae type k. If Z k is 1, B k can be any 
non-negative value; if it is 0, B k will be 0. This is modelled 
by adding three linear constraints for each algae type k. 


KL + KD + KW< K™ x + KM (1 - Z k ) (12.103) 

KL + KD + KW> Kf n (Z k ) (12.104) 

B k < BMZ k (12.105) 

where KM and BM are any large numbers no less than the 
largest possible value of the total extinction or biomass 
concentration, respectively. Since the objective of maxi- 
mizing the sum of all B k together with Equation 12.105 is 
to set each Z b value equal to 1 , it is only when the total 
extinction is outside of the extinction range K™ m to K™ x 
that the Z b value will be forced to 0. Equation 12. 105 then 
forces the corresponding B k to 0. 

5.9.4. Growth Limits 

For all algae types k, the maximum potential biomass 
concentration, B™ x (g dry/m 3 ), at the end of the time 
interval At (days) depends on the initial biomass concen- 
tration, B k , (g dry/m 3 ); the maximum gross production 
rate, Pg™ x (day -1 ); the respiration rate constant, 
R k (day -1 ); and a time and depth averaged production 
efficiency factor, E k . Mortality is not included in this 
computation. Using the net production rate constant, 
Pn k (= Pgjf ax E k — R k ) (day -1 ), for each algae type k: 

Bjf ax = b£ exp{Pn k Atj (12.106) 

This condition is taken into account by the optimization 
algorithm. For types with an initial biomass of zero, a 
small base level is used instead in order to allow 
previously absent species to start growing. 

5.9.5. Mortality Limits 

As in the case of growth, the mortality of each algae 
species is also constrained to prevent a complete removal 
within a single time step. The minimum biomass value of 
a species is obtained when there is no production, but 
only mortality. The minimum biomass, B™ (g dry/m 3 ), of 
type k at the end of time interval At depends on the 
initial biomass, B° (g dry/m 3 ), of type k and the specific 
mortality rate constant, M k (day -1 ), of type k. 

Bf n = B° k Z k exp{— M k Atj (12.107) 

These maximum and minimum values are computed for 
each individual algae type. However, the model sums each 
of these maximum and minimum values over all subtypes 
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within each species and applies it to the total biomass of 
the species. This way the maximum possible mortality 
cannot be exceeded, but transitions between limit types 
remain possible. 

X X X B r v species j 

k of species j k of species j k of species j 

(12.108) 

As mortality is computed according to a negative 
exponential function, the minimum biomass level is 
always positive; in other words, a species can never 
disappear completely. For computational purposes, the 
minimum biomass of a species is set to zero if it drops 
below some small threshold base level. It is conceivable 
that the amount of biomass, which should minimally be 
maintained according to Equation 12.108 exceeds the 
amount permitted according to the available amount of 
light energy (12.100). If this happens, Condition (12.100) 
is dropped. 

5.9.6. Oxygen-Related Processes 

The oxygen concentration in the water column depends 
in part on the biochemical and physical processes that 
either produce or consume oxygen. For the algae 
the model includes the production of oxygen (primary 
production), and its consumption through respiration. 
Oxygen is also consumed by the mineralization of 
detritus and other organic material (in the water column 
and bottom sediment), nitrification of ammonia to nitrate, 
and exchange of oxygen with the atmosphere (e.g. 
re-aeration). 

The mineralization of carbon detritus in the water 
column and bottom consumes oxygen at a molar ratio of 
1:1, equivalent to a ratio of 32/12 g 0 2 to 1 g C. The net 
growth of algae produces oxygen at a molar ratio of 1:1, 
equivalent to a ratio of 32/12 g 0 2 to 1 g C. The miner- 
alization of organic carbon in waste and carbon detritus 
in the water column and bottom consumes oxygen 
at a molar ratio of 1:1, equivalent to a ratio of 32/12 g 
0 2 to 1 g C. The nitrification reaction consumes oxygen 
at a molar ratio of 2:1, equivalent to a ratio of 64/14 g 
0 2 to 1 g N. For all algae, oxygen is produced during 
photosynthesis. The net growth of algae produces 
oxygen at a molar ratio of 1:1, equivalent to a ratio of 
32/12 g 0 2 to 1 g C. 


Bottom Oxygen 

The mineralization of organic material in the bottom 
sediment consumes oxygen, which must be supplied 
from the water column. The process consumes oxygen at 
a molar ratio (oxygen to carbon) of 1:1, equivalent to a 
ratio of 32/12 g 0 2 to 1 g C. 

Daily Oxygen Cycle 

Because oxygen is produced by the photosynthesis of 
algae during the daylight hours, there is a natural varia- 
tion of oxygen concentrations over the twenty-four-hour 
day-night cycle. Typically, oxygen concentrations are 
lowest in the early morning as oxygen is consumed 
during the night through the processes of algae respira- 
tion and organic material mineralization. During the 
daylight hours, oxygen is produced, and the highest 
values (often supersaturated) are typically found in the 
late afternoon. When biomasses are high, these variations 
may be large enough to cause low oxygen conditions 
during the night or in the early morning. 

In the traditional BLOOM calculation, the water 
quality processes are all calculated for a daily averaged 
situation. This is reflected by the choice of the input 
parameters for the light model: the daily averaged solar 
radiation and the day length. Reducing the time step would 
be the most straightforward way to include diurnal 
variations. The drawback, however, is a considerable 
increase in computation time. Thus an alternative approach 
has been adopted in the model. The total daily rate of 
primary production is computed first. Next, this produc- 
tion is distributed over the day. The model takes into 
account the day length, and oxygen production begins in 
the first daylight hour. Oxygen production increases during 
the morning, levels off at a (user-defined) maximum value 
for a period in the middle of the day, and decreases during 
the afternoon. There is no oxygen production during the 
night-time. The hourly oxygen production is combined 
with the daily averaged oxygen consumption processes 
and the re-aeration to produce an hourly value of oxygen 
concentration in the water. 

Maintenance Respiration 

Respiration in algae is a process in which organic carbon 
is oxidized, using oxygen to produce energy. The process 
occurs throughout the day and results in oxygen 
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consumption at a rate of 32 grams of oxygen per 
12 grams of carbon: 32/12 g 0 2 to 1 g C. Total respiration 
is divided into growth and maintenance respiration. 
Maintenance respiration is a first-order temperature 
dependent process. Species growth respiration is 
accounted for when calculating the net primary produc- 
tion. Maintenance respiration, and thus the amount of 
oxygen consumed, is governed by the temperature- 
dependent respiration rate constant, kjf s 6^ es(T-20) (day -1 ), 
the temperature, T (°C), of the water (higher respiration 
at higher temperatures), and the concentration of algae 
biomass, AlgC k . Each algae type can have a different 
respiration rate. Hence, the rate (gC 2 /mVday) at which 
carbon in algae is oxidized is 

dAlgC k /dt = k T k es % e5(T - 20) AlgC k (12.109) 

6. Simulation Methods 

Most of those who will be using water quality models will 
be using simulation models that are commonly available 
from governmental agencies (e.g. USEPA), universities, or 
private consulting and research institutions such as 
the Danish Hydraulics Institute, Wallingford software or 
WL | Delft Hydraulics (Ambrose et al. , 1996; Brown and 
Barnwell, 1987; Cerco and Cole, 1995; DeMarchi et al, 
1999; Ivanov et al., 1996; Reichert, 1994; USEPA, 2001; 
WL | Delft Hydraulics, 2003). 

These simulation models are typically based on 
numerical methods that incorporate a combination of 
plug flow and continuously stirred reactor approaches 
to pollutant transport. Users must divide streams, rivers, 
and lakes and reservoirs into a series of well-mixed 
segments or volume elements. A hydrological or hydro- 
dynamic model calculates the flow of water between all 
of these. In each simulation time step, plug flow enters 
these segments or volume elements from upstream 
segments or elements. Flow also exits from them to 
downstream segments or elements. During this time 
the constituents can decay or grow, as appropriate, 
depending on the conditions in those segments or 
volume elements. At the end of each time-step, the 
volumes and their constituents within each segment or 
element are fully mixed. The length of each segment 
or the volume in each element reflects the extent of 
dispersion in the system. 


6.1. Numerical Accuracy 

Water quality simulation models based on physical, 
biological and chemical processes typically include time 
rate of change terms such as dC/dt. While it is possible 
to solve some of these differential equations analytically, 
most water quality simulation models use numerical 
methods. The purpose of this section is not to explain 
how this can be done, but rather to point to some of the 
restrictions placed on the modeller because of these 
numerical methods. 

Consider first the relationship between the stream, 
river or lake segments and the duration of time steps, At. 
The basic first-order decay flux, dC/dt (g/m 3 /day), for a 
constituent concentration, C, that is dependent on a rate 
constant, k (day -1 ), is: 

dC/dt = -kC (12.110) 

The finite difference approximation of this equation can 
be written 

C(£ + At) - C(t) = - C(t)kAt (12.111) 

or 

C(t + At) = C(t)(l - kAt) (12.112) 

This equation can be used to illustrate the restriction 
placed on the term kAt. That term cannot exceed a value 
of 1 or else C(t+ At) will be negative. 

Figure 12.25 is a plot of various values of 
C(t + A t)/C(t) versus kAt. This plot is compared with 
the analytical solution resulting from the integration of 
Equation 12.110, namely: 

C(t + At) = C(t) exp{— kAt} (12.113) 

Reducing the value of At will increase the accuracy of the 
numerical solution. Hence, for whatever value of At, it 
can be divided by a positive integer n to become l/n th of 
its original value. In this case the predicted concentration 
C(t 4- At) will equal 

C(£ + At) = C(t)(l - kAt/nT (12.114) 

For example if kAt = 1, and n = 2, the final concentra- 
tion ratio will equal 

C(t + A t)/C(t) = (1 - 1/2) 2 = 0.25 (12.115) 

Compare this 0.25 to 0.37, the exact solution, and to 
0.0, the approximate solution when n is 1. Having n = 2 
brings a big improvement. If n = 3, the concentration 
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Figure 12.25. Plot of numerical approximation (red line) 
based on Equation 12.1 12 compared to the true analytical 
(blue line) value obtained from Equation 12.113. 


ratio will be 0.30, an even greater improvement compared 
to 0. However, no matter what value of n is selected, 
the predicted concentration will always be less than the 
actual value based on Equation 12.113, and hence the 
error is cumulative. Whenever At > n/fe, the predicted 
concentrations will alternate between positive and 
negative values, either diverging, converging or just 
repeating the cycle, depending on how much At exceeds 
n/k. In any event, the predicted concentrations are not 
very useful. 

Letting m = — n/kAt , Equation 12.114 can be 
written as 

C{t + At) = C(tXl + (12.116) 

As n approaches infinity, so does the variable m, and 
hence the expression (1 + l/m) m becomes the natural 
logarithm base e = 2.718282. Thus, as n approaches 
infinity, Equation 12.114 becomes Equation 12.113, the 
exact solution to Equation 12.110. 

6.2. Traditional Approach 

Most water quality simulation models simulate quality 
over a consecutive series of discrete time periods. Time is 
divided into discrete intervals and the flows are assumed 
constant within each of those time period intervals. Each 
water body is divided into segments or volume elements, 
and these are considered to be in steady-state conditions 
within each simulation time period. Advection or plug 
flow (i.e. no mixing or dispersion) is assumed during each 
time period. At the end of each period mixing occurs 


within each segment or volume element to obtain the 
concentrations in the segment or volume element at 
the beginning of the next time step. 

This method is illustrated in Figure 12.26. The indices 
i — 1 , i and i+1 refer to stream or river reach segments. 
The indices t and t+1 refer to two successive time peri- 
ods, respectively. At the beginning of time period t, each 
segment is completely mixed. During the time interval At 
of period t, the water quality model predicts the concen- 
trations, assuming plug flow in the direction of flow from 
segment i toward segment i + 1 . The time interval At is 
such that the flow from any segment i does not pass 
through any following segment i + 1. Hence, at the end of 
each time period each segment has some of the original 
water that was there at the beginning of the period, 
and its end-of-period concentrations of constituents, plus 
some of the immediately upstream segment’s water and 
its end-of-period concentrations of constituents. These 
two volumes of water and their respective constituent 
concentrations are then mixed to achieve a constant 
concentration within the entire segment. This is done for 
all segments in each time step. Included in this plug flow 
and then mixing process are the inputs to the reach from 
point and non-point sources of constituents. 

In Figure 12.26 a mass of waste enters reach i at a rate 
of W{. The volume in each reach segment is denoted by V, 
and the flows from one segment to the next are denoted 
by Q. The drawing on the left represents a portion of a 
stream or river divided into well-mixed segments. During 
a period £, waste constituents enter segment i from the 
immediate upstream reach i— 1 and from the point waste 
source. In this illustration, the mass of each of these 
wastes is assumed to decay during each time period, 
independent of other wastes in the water. Depending on 
the types of waste, the decay (or even growth) processes 
that take place may be more complex than those assumed 
in this illustration. At the end of each time period these 
altered wastes are mixed together to create an average 
concentration for the entire reach segment. This illustra- 
tion applies for each reach segment i and for each time 
period t. 

The length, Ax b of each completely mixed segment or 
volume element depends on the extent of dispersion. 
Reducing the length of each reach segment or size of each 
volume element reduces the dispersion within the entire 
stream or river. Reducing segment lengths, together with 
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Figure 12.26. Water quality 
modelling approach showing a 
water system schematized into 
computational cells. 



increasing flow velocities, also reduces the allowable 
duration of each time period t. The duration of each 
simulation time step At must be such that flow from any 
segment or element enters only the adjacent down- 
stream segment or element during that time step. Stated 
formally, the restriction is: 

At < T { (12.117) 

where T { is the residence time in reach segment or volume 
element i. For a 1-dimensional stream or river system 
consisting of a series of segments i of length Ax b cross 
section area A { and average flow Q it , the restriction is: 

At < min{AX|A/Q it ; Vi, tj (12.118) 

If time steps are chosen that violate this condition, then 
numerical solutions will be in error. The restriction 
defined by Equation 12.118 is often termed the 
‘courant condition’. It limits the maximum time-step 
value. Since the flows being simulated are not always 
known, this leads to the selection of very small time 
steps, especially in water bodies having very little 
dispersion. While smaller simulation time steps 
increase the accuracy of the model output, they also 


increase the computational times. Thus the balance 
between computational speed and numerical accuracy 
restricts the model efficiency in the traditional 
approach to simulating water quality. 

6.3. Backtracking Approach 

An alternative Lagrangian or backtracking approach to 
water quality simulation eliminates the need to consider 
the simulation time-step duration restriction indicated by 
Equation 12.118 (Manson and Wallis, 2000; Yin, 2002). 
The backtracking approach permits any simulation time- 
step duration to be used along with any segmenting 
scheme. Unlike the traditional approach, water can travel 
through any number of successive segments or volume 
elements in each simulation time step. 

This approach differs from the traditional one in that, 
instead of following the water in a segment or volume 
element downstream, the system tracks back upstream to 
find the source concentrations of the contaminants at 
time t that will be in the control volume or segment i + 1 
at the beginning of time period t + 1. 
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The backtracking process works from upstream to 
downstream. It starts from the segment of interest, i, and 
finds all the upstream sources of contaminants that 
flow into segment i during time period t. The contami- 
nants could come from segments in the same river reach 
or storage site, or from upstream river reaches or storage 
volume segments. They could also come from incremen- 
tal flows into upstream segments. Flows between the 
source site and the segment i+1 transport the contami- 
nants from their source sites to segment i during the time 
interval At, as shown in Figure 12.27. 

The simulation process for each segment and for 
each time period involves three steps. To compute the 
concentration of each constituent in segment i at the end 
of time period t, as shown in Figure 12.27, the approach 
first backtracks upstream to locate all the contaminant 
particles at the beginning of period t that will be in the 
segment i at the end of period t. This is achieved by 
finding the most upstream and downstream positions of 
all reach intervals that will be at the corresponding 
boundaries of segment i at the end of time period t. To do 


this requires computing the velocities through each of 
the intermediate segments or volume elements. Second, 
the changes in the amounts of the modelled quality 
constituents, such as temperature, organics, nutrients and 
toxics, are calculated, assuming plug flow during the 
time interval, At, and using the appropriate differential 
equations and numerical methods for solving them. 
Finally, all the multiple incoming blocks of water 
with their end-of-period constituent concentrations are 
completely mixed in the segment i to obtain initial 
concentrations in that segment for the next time step, 
t+1. This is done for each segment i in each time period 
t, proceeding in the downstream direction. 

If no dispersion is assumed, the backtracking process 
can be simplified to consider only the end points of each 
reach. Backtracking can take place to each end-of-reach 
location whose time of travel to the point of interest is just 
equal to or greater than At. Then, using interpolation 
between end-of-period constituent concentrations at 
those upstream sites, plus all loadings between those sites 
and the downstream site of interest, the constituent 



Figure 12.27. The backtracking 
approach for computing the 
concentrations of constituents in 
each reach segment or volume 
element / during time-step 
duration of At. 
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concentrations at the end of the time period t at the 
downstream ends of each reach can be computed. This 
process, like the one involving fully mixed reach 
segments, must take into account the possibility of 
multiple paths from each pollutant source to the site of 
interest, and the different values of rate constants, 
temperatures and other water quality parameters in each 
reach along those paths. 

Figure 12.27 illustrates an example of backtracking 
involving simple first-order decay processes. Assume 
contaminants that end up in reach segment i at the 
beginning of period t + 1 come from ] sources with initial 
concentrations C\, C\ , C 3 , ..., Cj at the beginning of time 
period t. Decay of mass from each source j during time 
At in each segment or volume element is determined by 
the following differential equation: 

dCj/df = - fe$ (T-20) q (12.119) 

The decay rate constant kp temperature correction coeffi- 
cient Oj and water temperature T are all temporally and 
spatially varied variables. Their values depend on the 
particular river reaches and storage volume sites through 
which the water travels during the period t from sites j to 
segment i. 

Integrating Equation 12.119 yields: 

C] +1 = C) exp{- (12.120) 

Since At is the time it takes water having an initial 
concentration Cj to travel to reach i, the values Cj +1 can be 
denoted as Qj +1 : 

cl +1 = Cj exp{— C# ! '“ 20 : A(} (12.121) 

In Equation 12.121 the values of the parameters are the 
appropriate ones for the stream or river between the 
source segments j and the destination segment i. These 
concentrations times their respective volumes, Vj, can 
then be mixed together to define the initial concentra- 
tion, C[ + 1 , in segment i at the beginning of the next 
time period t + 1 . 

6.4. Model Uncertainty 

There are two significant sources of uncertainty in water 
quality management models. One stems from incomplete 
knowledge or lack of sufficient data to estimate the 
probabilities of various events that might happen. 


Sometimes it is difficult to even identify possible future 
events. This type of uncertainty - sometimes called 
epistemic (Stewart, 2000) - stems from our incomplete 
conceptual understanding of the systems under study, by 
models that are necessarily simplified representations 
of the complexity of the natural and socio-economic 
systems, as well as by limited data for testing hypotheses 
and/or simulating the systems. 

Limited conceptual understanding leads to parameter 
uncertainty. For example, there is an ongoing debate 
about the parameters that can best represent the fate 
and transfer of pollutants through watersheds and water 
bodies. Arguably, more complete data and more work on 
model development can reduce this uncertainty. Thus, a 
goal of water quality management should be to increase 
the availability of data, improve their reliabilities and 
advance our modelling capabilities. 

However, even if it were possible to eliminate knowl- 
edge uncertainty, complete certainty in support of water 
quality management decisions will probably never be 
achieved until we can predict the variability of natural 
processes. This type of uncertainty arises in systems 
characterized by randomness. Assuming past observa- 
tions are indicative of what might happen in the future 
and with the same frequency - in other words, assuming 
stationary stochastic processes - we can estimate from 
these past observations the possible future events or 
outcomes that could occur and their probabilities. Even if 
we think we can estimate how likely any possible type of 
event may be in the future, we cannot predict precisely 
when or to what extent that event will occur. 

For ecosystems, we cannot be certain we know even 
what events may occur in the future, let alone their 
probabilities. Ecosystems are open systems in which it is 
not possible to know in advance what all the possible 
biological outcomes will be. Surprises are not only 
possible, but likely; hence, neither type of uncertainty 
- knowledge uncertainty nor unpredictable variability or 
randomness - can be eliminated. 

Thus, uncertainty is a reality of water quantity and 
quality management. This must be recognized when 
considering the results of water quality management 
models that relate to actions taken to meet the desired 
water quality criteria and designated uses of water bodies. 
Chapter 9 suggests some ways of characterizing this 
uncertainty. 
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7. Conclusions: Implementing a 
Water Quality Management Policy 

This chapter has provided only a brief introduction to 
some of the relationships contained in water quality 
models. As can be said for the other chapters as well, it 
summarizes a subject on which entire texts, and very 
good ones, have been written (see, for example, Chapra, 
1997; McCutcheon, 1989; Orlob, 1983; Schnoor, 1996; 
Thomann and Mueller, 1987). Water quality modelling 
and management demand skill and data. Skill comes with 
experience. Sufficient expertise will not be gained by 
working just with the material introduced in this chapter. 
It serves only as an introduction to surface water quality 
models, their assumptions and their limitations. 

If accompanied by field data and uncertainty analysis, 
many existing models can be used to assist those respon- 
sible for developing water quality management plans in 
an adaptive implementation or management framework. 
Adaptive implementation or management will allow 
for both model and data improvements over time. Such 
approaches strive toward achieving water quality stan- 
dards while relying on monitoring and experimentation 
to reduce uncertainty. This is often the only way one can 
proceed, given the complexity of the real world compared 
to the predictive models and the data and time usually 
available at the time a water quality analysis is needed. 
Starting with simple analyses and iteratively expanding 
data collection and modelling as the need arises is a 
reasonable approach. 

An adaptive management process begins with initial 
actions that have reasonable chances of succeeding. 
Future actions must be based on continued monitoring 
of the water body to determine how it responds to the 
actions taken. Plans for future regulation and public 
spending should be subject to revision as stakeholders 
learn more about how the system responds to actions 
taken. Monitoring is an essential aspect of adaptive water 
quality management and modelling (see Appendix B). 

Regardless of what immediate actions are taken, there 
may not be an immediate measurable response. For 
example, there may be significant lags between the time 
when actions are taken to reduce nutrient loads and 
the resulting changes in nutrient concentrations. This is 
especially likely if nutrients from past activities are tightly 


bound to sediments or if nutrient-contaminated ground- 
water has a long residence time before its release to 
surface water. For many reasons, lags between actions 
taken and responses must be expected. Water bodies 
should be monitored to establish whether the ‘trajectories’ 
of the measured water quality criteria point toward 
attainment of the designated use. 

Wasteload allocations will inevitably be required if 
quality standards are not being met. These involve costs. 
Different allocations will have different total costs and 
different distributions of those costs; hence they will have 
different perceived levels of fairness. A minimum-cost 
policy may result in a cost distribution that places most of 
the burden on just some of the stakeholders. But until 
such a policy is identified, one will not know this. An 
alternative may be to reduce loads from all sources by the 
same proportion. Such a policy has prevailed in the 
United States over the past several decades. Even though 
not very cost-effective from the point of view of water 
quality management, the ease of administration and the 
fulfilment of other objectives must have made such a 
policy politically acceptable, even though expensive. 
However, these types of wasteload allocations policies will 
not in themselves be sufficient for many of the ecosystem 
restoration efforts that are increasingly being made. 
Restoration activities are motivated in part by the services 
ecosystems provide for water quality management. 

Our abilities to include ecosystem components within 
water quantity and quality management models are at a 
fairly elementary level. Given the uncertainty, especially 
with respect to the prediction of how ecosystems will 
respond to water management actions, together with the 
need to take actions now, long before we can improve 
these capabilities, the popular call is for adaptive manage- 
ment. The trial and error aspects of adaptive management 
based on monitoring and imperfect models may not satisfy 
those who seek more definitive direction from water qual- 
ity analysts and their predictive models. Stakeholders and 
responsible agencies seeking assurances that the actions 
taken will always work as predicted may be disappointed. 
Even the best predictive capabilities of science cannot 
ensure that an action that will lead to the attainment of 
designated uses will be initially identified. Adaptive 
management is the only reasonable option in most cases 
for allowing water quality management programmes to 
move forward in the face of considerable uncertainties. 
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